REPORT  DOCUMENTATION  PAGE 


Form  Approved 
'  0MB  No.  074-0188 


Pitilic  reporting  liadon  for  ttwcolleao^inlbmiationteesliinajWtqavwagelhy  per  rMpcny,.»Tclu<lin9ljie^e  for  wvi^ 

needed  ari  complying  arKl  reviewing  this  colecl^  SerKi  comments  resting  this  biiden  estimate  or  any  other  asped  of  mtt<»»ectionrf  rtofirabon,hc^ 

bunten  to  VVbshsS^HeadnSrto^Sefvioes,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arfmgton,  VA  22202-4302,  and  to  the  Office  of  Managernent  and 
Budget  Paperwork  Redix:tion  Protect  (0704-0188).  Vifashingtoa  DC  20503  - 

1 .  AGENCY  USE  ONLY  (Leave  Wank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

April  1998  Final  Technical  Report  >  1/1/95  ^  9/30/97  _ 


4.  TITLE  AND  SUBTITLE  5.  FUNDING  NUMBERS 

LARGE  EDDY  SIMULATION  OF  SUPERSONIC  INLET  FLOWS  F-49620-95-1-0108 


6.  AUTHOR(S) 

PROF.  PARVIZ  MOBSr 
PROF.  SANJIVA  K.  LELE 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  AODRESS(ES) 

STANFORD  UNIVERSITY 
MECHANICAL  ENGINEERING 
FLOW  PHYSICS  &  COMPUTATION 
DIVISION 

STANFORD,  CA  94305-3030  _ 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


AFOSR 

DR.  LEONIDAS  SAKELL 
no  DUNCAN  AVE.,  STE.  B115 
BOLLING  AFB,  DC  20332-0001 

11.  SUPPLEMENTARY  NOTES 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

2DJA423 


10.  SPONSORING  /  MONITORING 
AGENCY  REPORT  NUMBER 


19980514  104 


12a.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

APPROVED  FOR  PUBLIC  RELEASE.  DISTRIBUTION  IS  UNLIMITED 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Max/fflum  200  Words;  .  •  •  ^  r,  m.- 

The  interaction  of  a  shock  wave  with  a  turbulent  boundary  layer  is  a  central  problem  m  supo'somc  mlet  flows. 
work  uses  numerical  and  analytical  techniques  to  study  shodotobulence  interaction  in  order  to  identify  and  explain 
factors  important  in  shock/boundary  layo*  interaction.  Direct  numerical  simulation  of  a  normal  shock  wave  with  an 
isotropic  turbulent  field  of  vorticity  and  entropy  fluctuations  showed  that  negative  upstream  correlation  between  tire 
vorticity  and  entropy  fluctuations  enhances  the  turbulence  across  the  shock.  Positive  upstream  correlation  has  a 
suppressing  effect.  A  new  numerical  method  providing  excellent  high  wavenumber  resolution  while  reducing  the 
computational  cost  was  developed.  A  model  with  no  ^justable  constants  was  developed  to  stady  the  vortex 
breakdown  resulting  from  the  interaction  of  canard  or  forbody  vortices  with  the  shock  waves  in  a  supersonic  inlet 
flow.  Very  good  agreement  with  both  experiment  and  computation  was  obtained.  A  numerical  method  to  compute 
shock/turbulence  interaction  using  a  conservative  form  of  Ae  Large  Eddy  Simulation  (LES)  equations  has  been 
developed  and  validated.  LES  of  the  interaction  of  isotropic  turbulence  with  a  normal  shock  was  performed  and 
comparisons  with  direct  numerical  simulation  (DNS)  results  were  favorable.  A  new  Fortran  90  code  has  been 
developed  for  the  computation  of  shock/turbulence  interaction.  The  code  is  an  improved  version  of  codes  used 
previously  in  shock/turbulence  interaction  simulations.  _ 


14.  SUBJECT  TERMS 

LARGE  EDDY  SIMULATION,  SUPERSONIC  INLET  FLOWS,  TURBULENCE, 
SHOCK  TURBULENCE 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

UNCLASSIFIED 


NSN  754(M)1-280-5500 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

UNCLASSIFIED 


19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

UNCLASSIFIED 


15.  NUMBER  OF  PAGES 
90 


16.  PRICE  CODE 


20.  LIMITATION  OF  ABSTRACT 


Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  Z39-18 
298-102 


Final  Technical  Report  for  the  period  Janu^  1, 1995  -  September  30, 1997  for  the  grant 

entitled: 


RJSEASCH  (/.FSC. 

’"©viewed  and  ic 


c 

D 

i;c>r-c.c,  . 

sum  Progrm  Msndci^y 


Large  Eddy  Simulation  of  Supersoiiit  Inlet  Flows 


Parviz  Moin  and  Sanjiva  K.  Lele 

Stanford  University 

Mechanical  Engineering,  Flow  Physics  &  Computation  Division 
Stanford,  CA  94305-3030 


Prepared  with  the  support  of  the 
Air  Force  Office  of  Scientific  Research 
under  AFOSR  Grant:  F49620-95-1-0108 


OVERVIEW 


This  report  summarizes  our  study  of  the  interaction  of  a  shock  wave  with  turbulent 
flows  for  the  period  1995-1997.  One  post-doctoral  fellow  and  one  graduate  student  were 
supported  by  this  program:  Dr.  Krishnan  Mahesh  obtained  his  degree  in  1996  (his  thesis 
was  published  as  Report  TF-69,  a  copy  was  sent  to  Dr.  Sakell)  and  Mr.  Albert  Honein  is 
continuing  his  doctoral  studies  in  the  Department  of  Mechanical  Engineering  at  Stanford. 
The  sahent  achievements  of  this  program  are  presented  below.  Details  are  provided  in  the 
attached  enclosures  of  relevant  reports. 

A  combination  of  linear  analysis  and  direct  simulation  was  used  by  Mahesh,  Moin 
and  Lele(^“^^  to  study  the  interaction  of  turbulent  shear  flow  with  a  normal  shock  wave. 
Subsequent  to  graduation,  Dr.  Mahesh  continued  his  work  as  a  post-doctoral  fellow  at  the 
Center  for  Turbulence  Research.  His  principal  accomplishments  are  as  follows: 

1.  Direct  numerical  simulation  was  used  to  study  the  interaction  of  a  normal  shock 
wave  with  an  isotropic  turbulent  field  of  vorticity  and  entropy  fluctuations.  The 
upstream  correlation  between  the  vorticity  and  entropy  fluctuations  was  shown  to 
strongly  influence  the  evolution  of  the  turbulence  across  the  shock.  The  validity  of 
Morkovin’s  hypothesis  behind  a  shock  was  also  examined;  shock-front  oscillations  were 
found  to  invalidate  the  part  of  the  hypothesis  relating  Urms  and  Trms-  For  further 
details,  refer  to  Part  4  of  the  report. 

2.  A  new  numerical  method  providing  excellent  high  wavenumber  resolution  while  reduc¬ 
ing  the  computational  cost  was  developed^^).  This  method  is  now  being  incorporated 
in  the  shock/boxmdary  layer  codes.  Refer  to  Part  3  for  details. 

3.  The  vortex  breakdown  resulting  from  the  interaction  of  canard  or  forbody  vortices 
with  the  shock  waves  in  a  supersonic  inlet  flow  was  studied.  A  model  with  no  ad¬ 
justable  constants  was  developed^®)  and  very  good  agreement  with  both  experiment 
and  computation  was  obtained.  Refer  to  part  2  for  details. 

The  project  was  extended  by  Mr.  Honein  to  develop  the  technology  of  large-eddy 
simulation  (LES)  for  turbulent  flows  with  shock  waves.  His  principal  achievements  are 
summarized  below: 

1.  A  numerical  method  to  compute  shock/turbulence  interaction  using  LES  has  been 
developed.  The  important  feature  of  this  method  is  that  the  energy  equation  is  written 


in  a  conservative  form,  which  is  essential  in  computing  flows  with  shock  waves.  The 
method  has  been  validated  by  computing  spatially  decaying  compressible  turbulence. 

2.  The  interaction  of  isotropic  turbulence  with  a  normal  shock  was  computed.  For  low 
values  of  Rex,  comparisons  with  direct  numerical  simulation  (DNS)  results  were  favor¬ 
able.  The  trends  found  with  DNS  at  low  Rex  were  also  observed  with  LES  at  higher 
Rex. 

3.  A  new  code  has  been  developed  for  the  computation  of  shock/turbulence  interaction. 
The  code  is  an  improved  version  of  codes  used  previously  in  shock/turbulence  inter¬ 
action  simulations.  It  is  currently  being  used  to  study  the  interaction  of  a  flat  plate 
boundary  layer  with  a  normal  shock  wave. 

This  report  is  organized  as  follows.  The  development  of  LES  for  shock/turbulence 
interaction  is  described  in  Part  1.  Part  2  deals  with  the  model  used  to  predict  the  onset 
of  shock-induced  vortex-breakdown.  The  new  numerical  method  is  described  in  Part  3. 
Part  4  concludes  the  report  with  the  influence  of  entropy  fluctuations  on  the  interaction 
of  turbulence  with  shock  waves. 
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PART  1 


Large  eddy  simulation  of 
shock  turbulence  interaction 


1.  Motivation  and  objectives 

The  presence  of  shock  waves  is  an  essential  characteristic  of  high  speed  flows.  Under¬ 
standing  the  mechanisms  of  turbulence  interacting  with  a  shock  wave  is  of  fundamental 
importance  in  predicting  the  interaction  of  turbulent  boundary  layers  with  shock  waves 
which  occur  in  many  engineering  apphcations.  The  high  Reynolds  number  in  practical 
flows  makes  direct  numerical  simulation  (DNS)  unfeasible.  Large  eddy  simulation  (LES), 
where  the  large  scales  of  turbulence  are  resolved  and  small  scale  effects  are  modeled,  is 
thus  required  for  detailed  analysis. 

The  goal  of  this  work  is  to  study  shock/wave  boundary  layer  interactions  using  LES. 
Toward  this  goal  a  simpler  problem  of  the  spatial  evolution  of  compressible  isotropic  tur¬ 
bulence  and  its  interaction  with  a  shock  wave  were  studied  first.  The  developed  numerical 
method  for  LES  was  then  incorporated  in  a  new  code  for  the  simulation  of  shock/boundary 
layer  interaction.  The  simulations  conducted  to  test  and  develop  the  numerical  method 
for  LES  are  described  below.  A  description  of  the  new  code  is  also  provided. 

2.  Accomplishments 

A  nonconservative  formulation  of  the  energy  equation  (internal  energy)  was  used  to 
perform  large  eddy  simulations  of  compressible  turbulence  in  Moin  et  al.  (1991).  A 
conservative  form  of  the  energy  equation  is  highly  desirable  for  computing  flows  containing 
shocks.  A  conservative  set  of  equations  for  LES  were  derived  from  the  nonconservative 
equations  derived  by  Moin  et  al.  (1991).  The  validation  of  the  formulation  for  spatially 
decaying  turbulence  is  described  in  §2.1.  Performance  of  the  formulation  was  compared 
with  a  DNS  of  isotropic  turbulence/shock  wave  interaction,  which  is  reported  in  §2.2.  The 
new  code  for  shock/turbulence  interaction  is  described  in  §2.3. 

2.1  Spatially  evolving  compressible  turbulence 

A  separate  temporal  simulation  is  advanced  in  time  until  realistic  turbulent  condi¬ 
tions  are  realized.  The  resulting  field  is  used  to  generate  inflow  turbulence  for  the  spatial 
calculation.  The  simulations  are  advanced  in  time  until  statistical  convergence  is  reached. 


Taylor’s  hypothesis  is  then  invoked  to  compare  the  spatial  simulation  to  the  corresponding 
temporal  one.  Comparisons  between  the  two  simulations  are  shown  in  Figures  1,  2,  and 
3.  The  computation  corresponds  to  a  microscale  Reynolds  number  B.e\  =  24.4  and  a  tur¬ 
bulent  Mach  number  Mt  =  0.39,  specified  at  the  infiow  of  the  spatial  simulation.  A  mesh 
with  33^  grid  points  was  used.  As  seen  in  Figures  1,  2,  and  3,  good  agreement  is  obtained. 

For  higher  values  of  Rex,  it  was  noted  that  a  build-up  of  aliasing  errors  and  insuf¬ 
ficient  model  dissipation  prevented  the  simulations  from  being  advanced  for  long  times. 
Increasing  the  streamwise  number  of  points,  and  using  the  skew  ssrmmetric  form  in  dis¬ 
cretizing  the  continuity  equation  helped  in  stabilizing  the  simulations.  Similar  agreement 
with  the  corresponding  temporal  simulation  was  found  except  for  the  evolution  of  the 
density  fluctuations. 

2.2  Isotropic  turbulence/shock  wave  interaction 

The  LES  formulation  was  next  applied  to  the  interaction  of  isotropic  turbulence  with 
a  shock  wave.  The  paraaneters  used  were  chosen  to  reproduce  a  DNS  of  isotropic  ttu-bu- 
lence/shock  wave  interaction  performed  by  Mahesh  et  al.  (1997).  The  mean  Mach  nmnber 
was  1.29  and  the  inflow  turbulence  had  a  Rex  of  19.1  and  a  turbulent  Mach  munber  M* 
of  0.14. 

A  summary  of  the  numerical  method  used  follows.  A  uniform  mesh  is  used  in  the  two 
homogeneous  directions  and  a  stretched  mesh  is  used  in  the  non-homogeneous  direction. 
In  computing  spatial  derivatives,  a  sixth  order  essentially  non-oscillatory  (ENO)  shock 
capturing  scheme  is  used  for  the  shock  region  while  a  sixth  order  Fade  scheme  is  used 
outside  of  this  region.  A  third  order  Rimge-Kutta  method  is  used  for  time  advancement. 
Finally,  a  sponge  layer  where  the  flow  is  forced  towards  the  laminar  solution  is  used  for 
the  outflow  boundary  condition. 

Resiilts  shown  in  Figures  4  and  5  indicate  that  shock/tinbulence  interaction  at  this 
Rex  is  well  predicted  by  LES,  with  moderate  savings  in  computation  time  and  in  memory 
usage.  Simulations  at  a  higher  inflow  Rex  of  49,  where  DNS  cannot  be  used,  were  also 
performed  (Figures  6  and  7).  It  is  observed  that  the  trends  found  at  low  Rex  are  also 
obtained  at  higher  Rca. 

2.3  New  code  for  shock/turbulence  interaction 

A  spectral,  finite-difference  hybrid  code  was  developed  and  is  currently  being  applied 
to  the  normal  shock  wave/flate  plate  boundary  layer  interaction  calculations.  The  nu¬ 
merical  methodology  used  in  our  previous  computations  of  shock/turbulence  interaction 
was  implemented.  In  addition,  a  spectral  Fourier  collocation  scheme  was  employed  in  the 


homogeneous  spanwise  direction.  This  code  is  written  in  Fortran  90  and  is  an  optimized 
version  of  our  previous  codes  used  in  the  shock/turbulence  calculations.  It  runs  foiu  times 
faster  and  requires  only  half  the  memory.  Finally,  it  is  highly  structured  and  will  be  easily 
ported  and  optimized  for  use  on  parallel  machines. 
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Figure  1.  Comparison  of  the  temporal  ( - )  and  spatial  (•  )  evolution  of 

the  turbulent  kinetic  energy.  Conditions  at  the  inflow  correspond  to  Rex=2AA  and 
Mt=0.39  {jt  is  the  turbulence  time  scale). 
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Figure  2.  Comparison  of  the  temporal  ( - )  and  spatial  (•  )  evolution  of  the 

density  fluctuations.  Conditions  at  the  inflow  correspond  to  Rex=2AA  and  Mt=0.39 
(tj  is  the  turbulence  time  scale). 
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Figure  3.  Comparison  of  the  temporal  ( - )  and  spatial  (•  )  evolution  of 

the  pressure  fluctuations.  Conditions  at  the  inflow  correspond  to  Re\=2AA  and 
Mi=0.39  (r*  is  the  turbulence  time  scale). 
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Figure  4.  Streamwise  evolution  of  the  turbulent  kinetic  energy  components  across 
•  a  Mach  1.29  shock  wave  with  inflow  Re^  =  19.1.  All  the  curves  are  normalized  by 

their  value  immediately  upstream  of  the  shock.  DNS  :  m  x  v'^;  LES  :  . u'^, 
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Figure  5.  Streamwise  evolution  of  vorticity  fluctuations  across  a  Mach  1.29 
shock  wave  vdth  inflow  Rex  =  19.1.  All  the  curves  are  normalized  by  their  value 
immediately  upstream  of  the  shock.  DNS  :  •  w^,  x  W2  =  LES  : .  w^, 


Figure  6.  Streamwise  evolution  of  the  turbulent  kinetic  energy  components  from 
LES  of  a  Mach  1.29  shock  wave  with  inflow  Rex  =  49.  All  ^  curves  ^  normalized 
by  their  value  immediately  upstream  of  the  shock. - 
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Figure  7.  Streamwise  evolution  of  vorticity  fluctuations  from  LES  of  a  Mach 
1.29  shock  wave  with  inflow  Rex  =  49.  All  the  curves  are  normalized  by  their  value 
immediately  upstream  of  the  shock.  - , - W2  =  • 


PART  2 


A  model  for  the  onset  of  breakdown  in  an  axisymmetric 
compressible  vortex 

Krishnan  Mahesh®^ 

Center  for  Turbulence  Research,  Stanford  University,  Stanford,  California  94305 
(Received  28  May  1996;  accepted  5  September  1996) 

A  simple  inviscid  model  to  predict  the  onset  of  breakdown  in  an  axisymmetric  vortex  is  proposed. 
Three  problems  are  considered:  the  shock-induced  breakdown  of  a  compressible  vortex,  the 
breakdown  of  a  free  compressible  vortex,  and  the  breakdown  of  a  free  incompressible  vortex.  The 
same  physical  reasoning  is  used  in  all  three  problems  to  predict  the  onset  of  breakdown.  It  is 
hypothesized  that  breakdown  is  the  result  of  the  competing  effects  of  adverse  pressure  rise  and 
streamwise  momentum  flux  at  the  vortex  centerline.  Breakdown  is  assumed  to  occur  if  the  pressure 
rise  exceeds  the  axial  momentum  flux.  A  formula  with  no  adjustable  constants  is  derived  for  the 
critical  swirl  number  in  all  three  problems.  The  dependence  of  the  critical  swirl  number  on 
parameters  such  as  upstream  Mach  number,  excess/deficit  in  centerline  axial  velocity,  and  shock 
oblique  angle  is  explored.  The  predictions  for  the  onset  of  shock-induced  breakdown  and  free 
incompressible  breakdown  are  compared  to  experiment  and  computation,  and  good  agreement  is 
observed.  Finally,  a  new  breakdown  map  is  proposed.  It  is  suggested  that  the  adverse  pressure  rise 
at  the  vortex  axis  be  plotted  against  the  axial  momentum  flux  to  determine  the  onset  of  breakdown. 
The  proposed  map  allows  the  simultaneous  comparison  of  data  from  flows  ranging  from 
incompressible  breakdown  to  breakdown  induced  by  a  shock  wave.  ©  1996  American  Institute  of 
Physics.  [81070-6631(96)02112-5] 


I.  INTRODUCTION 

A  large  body  of  information  exists  (e.g.  see  the  review 
articles  by  Hall,^  Leibovich,^  Wedemeyer,^  Escudier,"^ 
Stuart,^  and  Delery®)  on  the  breakdown  of  incompressible 
streamwise  vortices.  Less  is  known  about  vortex  breakdown 
at  high  speeds.  An  interesting  example  of  supersonic  vortex 
breakdown  is  the  breakdown  induced  by  the  interaction  of 
the  vortex  with  a  shock  wave.  The  flow  in  supersonic  engine 
inlets  and  over  high-speed  delta  wings  constitute  technologi¬ 
cally  important  examples  of  this  phenomenon,  which  is 
termed  “shock-induced  vortex  breakdown.” 

Gustintsev  etalJ^  and  Zatoloka  etal^  appear  to  have 
conducted  the  earliest  investigations  into  shock-induced  vor¬ 
tex  breakdown.  Tlie  qualitative  similarity  of  the  flow  to  that 
of  a  separated  boundary  layer  was  noted  in  these  experi¬ 
ments.  Subsequently,  Horowitz,^  Delery  et  Metwally 
etal}^  and  Cattafesta  and  Settles^^  have  experimentally 
studied  vortex  breakdown  induced  by  a  normal  shock.  The 
interaction  between  streamwise  vortices  and  wedge-attached 
oblique  shock  waves  was  experimentally  investigated  by 
Kalkhoran  and  Sforza.^^ 

Horowitz^  and  Delery  et  al^^  were  the  first  to  quantita¬ 
tively  characterize  the  nature  of  the  breakdown.  Their  experi¬ 
ments  studied  normal  shocks  of  strength  equal  to  Mach  1.6, 
1.75, 2  and  2.28.  At  each  Mach  number,  they  varied  the  swirl 
in  the  incident  vortex  and  identified  a  critical  swirl  number 
above  which  the  vortex  would  break  down.  The  results  were 
plotted  on  a  “breakdown  map”  of  swirl  number  against 
Mach  number,  where  it  was  observed  that  the  critical  swirl 
number  decreased  as  the  Mach  number  of  the  shock  in¬ 
creased.  A  companion  numerical  study  using  the  Euler  equa- 
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tions  supported  the  experimentally  observed  trends.  The  ex¬ 
periments  by  Metwally  etal}^  and  Cattafesta  and  Settles^^ 
extended  the  range  of  available  data  to  Mach  4.  Based  on 
their  visualization  of  the  flow,  Metwally  et  a/.“  proposed  a 
qualitative  picture  of  the  flow-field  resulting  from  the  break¬ 
down  of  the  vortex. 

Rizzetta^"^  obtained  numerical  solutions  to  the  Reynolds 
averaged  Euler  and  Navier-Stokes  equations,  with  the  objec¬ 
tive  of  predicting  Kalkhoran  and  Sforza’s^^  experimental 
measurements  of  pressure  distribution  on  the  wedge.  The 
swirling  supersonic  flow  in  a  circular  duct  was  computed  by 
Kandil  etal}^'^^  who  provided  qualitative  flow-field  infor¬ 
mation  on  the  breakdown.  The  most  extensive  computations 
of  shock-induced  vortex  breakdown  are  the  recent  calcula¬ 
tions  by  Erlebacher  et  al}^  These  workers  studied  the  inter¬ 
action  between  a  streamwise  vortex  and  a  normal  shock 
wave  using  the  unsteady,  axisymmetric,  compressible 
Navier-Stokes  equations.  Mach  numbers  from  1.3  to  10 
were  computed.  In  the  same  spirit  as  Delery  et  a  critical 
swirl  number  was  numerically  identified  at  each  Mach  num¬ 
ber,  and  a  breakdown  map  of  swirl  number  against  Mach 
number  made.  The  trend  observed  by  Delery  etal}^  was 
seen  to  extend  to  Mach  10;  i.e.,  the  critical  swirl  number 
decreased  with  increasing  Mach  number.  Some  interesting 
features  of  the  flow  field  were  also  highlighted. 

The  only  attempt  to  quantitatively  predict  some  aspect  of 
shock-induced  breakdown  appears  to  have  been  made  by 
Cattafesta^*  who  equated  the  ratio  of  swirl  number  (down¬ 
stream  to  upstream)  across  the  shock  wave  to  the  velocity 
ratio  (upstream  to  downstream)  across  the  shock.  By  com¬ 
paring  to  experimental  data,  he  obtained  a  value  of  0.6  for 
the  swirl  number  behind  the  shock  wave.  More  recently,  Er¬ 
lebacher  etal}^  have  proposed  an  empirical  correlation  be¬ 
tween  the  critical  swirl  number  and  the  Mach  number  of  the 
shock  wave,  based  on  a  curve  fit  to  their  data. 
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FIG.  1.  Schematic  of  the  interaction  between  a  strcamwise  vortex  and  a 
normal  shock  wave. 


In  this  paper,  we  propose  a  model  to  predict  the  onset  of 
shock-induced  vortex  breakdown.  The  proposed  model  has 
no  adjustable  constants,  and  is  compared  to  both  experiment 
and  computation.  Also,  the  dependence  of  the  critical  swirl 
number  on  parameters  such  as  the  upstream  Mach  number, 
excess/deficit  in  centerline  axial  velocity,  and  shock  oblique 
angle  is  explored.  Two  other  problems  are  then  considered: 
the  breakdown  of  a  free  compressible  vortex,  and  free  in¬ 
compressible  vortex  breakdown.  The  same  breakdown  crite¬ 
rion  is  used  in  all  three  problems  to  predict  the  onset  of 
breakdown.  Finally,  a  new  breakdown  map  is  proposed,  that 
allows  the  simultaneous  comparison  of  data  from  flows  rang¬ 
ing  from  incompressible  breakdown  to  breakdown  induced 
by  a  shock  wave. 

This  paper  is  organized  as  follows.  A  description  of  the 
upstream  vortex  is  first  provided  in  Section  IIA.  This  is 
followed  in  Section  n  B  by  a  description  of  the  proposed 
breakdown  criterion  and  expressions  for  the  critical  swirl 
number.  Section  m  compares  the  model  predictions  to  com¬ 
putation  and  experiment.  The  influence  of  centerline  excess/ 
deficit  in  axial  velocity,  and  obliquity  of  the  shock  wave  is 
also  discussed.  The  onset  of  breakdown  in  a  free  compress¬ 
ible  vortex  is  discussed  in  Section  IV.  Incompressible  vortex 
breakdown  is  briefly  considered  in  Section  V.  A  new  break¬ 
down  map  is  then  proposed  in  Section  VI.  The  paper  is  con¬ 
cluded  with  a  brief  summary  in  Section  Vn. 

II.  STATEMENT  OF  PROBLEM 

Figure  1  shows  a  schematic  of  the  interaction  between  a 
stream^se  vortex  and  a  normal  shock  wave.  The  axial  flow 
is  from  left  to  right.  The  variables  x  and  r  are  used  to  denote 
the  axial  and  radial  coordinate  respectively.  The  axial  and 
swirl  components  of  velocity  are  denoted  by  U  and  re¬ 
spectively,  and  p,  p  and  T  represent  the  pressure,  density  and 
temperature.  The  subscripts  “oo”  and  “c”  correspond  to 
values  in  the  free-stream  and  the  centerline  of  the  vortex,  and 
the  states  upstream  and  downstream  of  the  shock  wave  are 
respectively  denoted  by  the  subscripts  “T’  and  “2”  (e.g., 
Poc2  denotes  the  free-stream  pressure  downstream  of  the 
shockwave). 

A  description  of  the  incident  vortex  is  first  provided  in 
Section  II  A.  This  is  followed  in  Section  II B  by  an  outline  of 
the  model. 


A.  The  upstream  vortex 

Studies  of  incompressible  vortex  breakdown  (e.g. 
Darmofal*^)  suggest  that  the  onset  of  breakdown  is  generally 
independent  of  viscosity  for  vortex  Reynolds  number  (based 
on  free-stream  axial  velocity  and  core  radius)  greater  than 
about  300.  As  a  result,  viscosity  is  neglected  in  this  paper. 
The  upstream  vortex  is  therefore  governed  by  the  axisym- 
metric,  compressible  Euler  equations.  It  is  readily  seen  that 
the  profiles. 


U=(/(r),  p=p(r),  p=p(r)  (1) 


trivially  satisfy  the  continuity,  axial  momentum  and  energy 
equations.  The  radial  momentum  equation. 


dp  _  pvl 
dr  r 


(2) 


remains  to  be  satisfied.  Experiments^®’^^  show  that  the  swirl 
profile  of  the  Burgers  vortex  is  a  good  fit  to  experimental 
data.  However,  the  Burgers  profile  makes  analytical  solution 
difficult.  As  a  result,  this  paper  uses  the  Rankine  vortex  as  an 
approximation  for  the  upstream  vortex.  Non- 
dimensionalizing  the  radial  coordinate  by  the  core  radius  (lo¬ 
cation  where  u  ^  is  maximum)  and  velocity  by  tiie  peak  value 
of  the  swirl  velocity  (denoted  by  u^),  the  swiri  velocity 
profile  of  the  upstream  vortex  is  given  by. 


Ve-n  r^l 


(3) 


where  the  tilde  is  used  to  denote  non-dimensional  variables. 

The  density  varies  with  radius  for  a  compressible  vortex. 
This  paper  considers  two  different  idealizations  of  the  ther¬ 
modynamic  field  in  the  upstream  vortex:  spatially  uniform 
stagnation  temperature  and  spatially  uniform  entropy.  The 
assumption  of  uniform  stagnation  temperature  is  prompted 
by  experimental  data.  Delery  et  note  fliat  the  total  tem¬ 
perature  in  the  upstream  vortex  in  their  experiments  is  ap¬ 
proximately  uniform.  Measurements  in  a  Mach  3  vortex  by 
Metwally  etaO^  and  Cattafesta  and  Settles*^  seem  to  sup¬ 
port  this  approximation.  Cattafesta  and  Settles’  data  (Fig.  7 
of  their  paper)  show  a  deficit  of  about  4%  of  free-stream  in 
total  temperature  at  the  centerline.  The  idealization  of  uni¬ 
form  entropy  is  prompted  by  past  theoretical  and  computa¬ 
tional  studies  on  compressible  vortices  (e.g.  Colonius 
et  ai^. 

Expressions  for  the  centerline  pressure  and  density  for 
the  uniform  stagnation  temperature  vortex  and  uniform  en¬ 
tropy  vortex  are  obtained  below.  Defining  the  non- 
dimensional  variables. 


the  radial  momentum  equation  becomes. 


dr 


y(rAf„)2p-^. 


(4) 


(5) 


Phv^Fluids^Jo!^JJ^2jDecember1996 


Krishnan  Mahesh  3339 


The  variable  y  denotes  the  ratio  of  specific  heats  and  is  taken 
as  1.4  in  this  paper.  F  is  the  swirl  number  of  the  vortex,  and 
is  defined  as  F = t> /f/« .  Af „  is  the  free-stream  Mach  num¬ 
ber,  defined  as  M^=  UJc^.  F M,  will  be  recognized  as  the 
swirl  Mach  number,  V0„,lc^. 

Uniform  entropy  vortex:  If  the  entropy  is  spatially  uni¬ 
form, 

p—^~  (6) 

Expressing  the  density  in  terms  of  the  pressure  in  the  radial 
momennim  equation  and  integrating  yields  the  following  ex¬ 
pressions  for  the  centerline  pressure  and  density; 

Pc=[i-(y-i)r"A/2 

Uniform  stagnation  temperature  vortex:  The  spatial  uni¬ 
formity  of  sugnation  temperature  requires  that 


r+ 


ul 

- ^ - 

2Cp  2C/ 


Delery  et  experiments  show  that  the  axial  velocity  in 
the  upstream  vortex  was  nearly  uniform;  i.e.,  £/« .  Cat- 

tafesta  and  Settles*^  on  the  other  hand,  observe  a  wake-like 
profile.  This  paper  assumes  uniform  axial  velocity  for  the 
uniform  stagnation  temperature  vortex.  This  yields  the  fol¬ 
lowing  expression  for  the  non-dimensional  temperature  in 
the  vortex: 


f=l-^(FM„)V/.  (9) 

The  equation  of  state  implies  that  p=pT.  Substituting  for 
p  and  r  in  the  radial  momentum  equation  and  integrating 
yields  the  following  expressions  for  the  density  and  pressure 
at  the  centerline  of  the  uniform  stagnation  temperature  vor¬ 
tex: 


2y/(y-l) 


-  _  -  (10) 

Tf  —  1,  Pc~Pc‘ 

B.  A  criterion  for  shock-induced  breakdown 

A  simple  criterion  for  breakdown  of  the  upstream  vortex 
is  first  proposed.  The  properties  of  the  upstremn  vortex  (Sec¬ 
tion  n  A)  are  then  us^  to  obtain  an  expression  for  the  criti¬ 
cal  swirl  number  above  which  the  vortex  would  break  down. 
The  breakdown  criterion  is  based  upon  an  approximation  to 
the  axial  momentum  equation  at  the  centerline  of  the  vortex. 
Note  that  as  a  result  of  axisymmetry,  the  radial  velocity  at 
the  centerline  would  be  zero.  When  combined  with  the  swirl 
velocity  being  zero  at  the  centerline,  this  suggests  that  the 
flow  near  the  vortex  centerline  would  largely  be  in  the 
streamwise  direction.  The  one-dimensional  momentum  equa¬ 
tions  may  therefore  be  used  to  model  the  flow  around  the 
vortex  centerline.  p+pU^  would  therefore  be  constant 
across  a  region  of  rapid  streamwise  variation. 

Consider  the  vortex  impinging  upon  the  shock  wave.  On 
account  of  the  rotation,  the  pressure  at  the  center  of  the  vor- 


Pc= 


1  - 


tex  is  less  than  the  free-stream  value;  i.e., Pc\<Pm\-  Pressure 
rises  across  a  shock  wave;  i.e.,  Poo2>Pot>|.  The  vortex  there¬ 
fore  experiences  an  adverse  streamwise  pressure  rise,  which 
may  be  quantified  by  the  pressure  difference,  Pco2~Pci-  The 
fluid  in  the  vortex  has  a  certain  inertia  in  Ae  streamwise 
direction,  which  may  be  quantified  by  the  streamwise  mo¬ 
mentum  flux,  PciUh.  Breakdown  is  assumed  to  occur  if  the 
axial  pressure  rise  exceeds  the  upstream  streamwise  momen¬ 
tum  flux,  thereby  stagnating  the  flow;  i.e.,  if 

P-2-Pci^Pcit/?i^Pcif/i,(l  +  |^)  (11) 

where  AI/  denotes  the  upstream  excess  in  axial  velocity  at 
the  centerline.  If  the  axial  veloci^  is  uniform,  then 
Af/=0.  The  threshold  for  breakdown  is  therefore  given  by 
the  relation, 

,  /  Al/\^ 

P«>2-Pci=PciUiiy  +  jj^j  .  (12) 

The  axial  velocity  is  assumed  to  be  uniform  through  most  of 
this  paper.  The  effect  of  non-uniform  axial  velocity  is  sepa¬ 
rately  discussed  in  Section  m  B.  Equation  (12)  may  be  re¬ 
written  in  non-dimensional  form  for  uniform  axial  velocity 
as. 


2 

Poo2  “Pc  1 = yPciAf  Ml- 


(13) 


We  have  already  obtained  expressions  for  pd  and  p^  in 
terms  of  F  and  The  Rankine-Hugoniot  equations  for  a 
normal  shock  express  pv,2  in  terms  of  the  upstream  Mach 
number,  M«,i.  Substituting  for  Pd,  Pd  3®<1  ^®2  “to  tl*o 
above  breakdown  criterion  will  therefore  yield  an  expression 
for  the  critical  swirl  number  Fcnt  “  terms  of  Mach  number  of 
the  shock  wave  for  a  vortex  with  uniform  axial  velocity.  This 
expression  is  derived  below. 

Uniform  stagnation  temperature  vortex:  For  a  uniform 
stagnation  temperature  vortex,  we  have  pd~Pc\-  Substitu¬ 
tion  into  the  criterion  for  breakdown  [Eq.  (13)]  yields. 


P«.2 

i  +  yM^,’ 


(14) 


where  p»2  is  given  by  the  Rankine-Hugoniot  equations  as, 


i+^(wi,-i). 


(15) 


Substituting  for  Pd  from  Eq.  (10)  and  ^2  front  Eq.  (15)  into 
Eq.  (14),  we  get. 


1- 


1 


2y/(y-I) 


1  +  yM 


(16) 


which  upon  rearrangement  yields  the  following  expression 
for  the  critical  swirl  number  as  a  function  of  the  Mach  num¬ 
ber  of  the  shock: 
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Uniform  entropy  v£>rtex;  Expressions  for  the  centerline 
density  and  temperature  for  a  uniform  entropy  vortex  are 
given  by  Eq.  (7).  Substitution  into  Eq.  (13)  yields  the  follow¬ 
ing  implicit  equation  for  the  critical  swirl  number  as  a  func¬ 
tion  of  the  Mach  number: 

1  -(r- 

(is) 

The  Newton— Raphson  method  was  used  to  solve  the  above 
equation  for  Tent  as  a  function  of  the  Mach  number  of  the 
shock  wave. 

III.  RESULTS:  SHOCK-INDUCED  VORTEX 
BREAKDOWN 

A.  Uniform  axial  velocity 

Results  for  the  critical  swirl  number  are  presented  for  the 
case  where  the  axial  velocity  is  uniform.  Figure  2  shows  the 
predicted  values  of  the  critical  swirl  number  as  a  function  of 
the  Mach  number  of  the  shock.  The  predicted  values  are 
compared  to  the  experimental  values  reported  by  Delery 
et  (the  data  were  obtained  from  Fig.  35  of  their  paper) 
for  Mach  numbers  of  1.75,  2  and  2.28.  Also  shown  are  re¬ 
sults  from  the  computations  by  Erlebacher  etal}^  (the  data 
were  obtained  from  Table  3  of  their  report).  Note  that  the 
computational  data  at  Mach  1.7  were  very  close  to  the  ex¬ 
perimental  value  at  Mach  1.75  (0.331  as  compared  to  0.33). 
This  made  tiie  experimental  data  hard  to  discern  when  both 
experimental  and  computational  results  were  plotted.  As  a 


result,  the  computational  value  at  Mach  1.7  is  not  plotted  in 
Fig.  2. 

The  predicted  values  are  seen  to  be  in  good  agreement 
with  both  experiment  and  computation.  The  critical  swirl 
number  is  predicted  to  decrease  with  increasing  Mach  num¬ 
ber  as  observed.  According  to  the  proposed  criterion  [Eqs. 
(7),  (10)  and  (13)],  this  decrease  in  Fcrit  is  due  to  a  combi¬ 
nation  of  two  factors:  increase  in  the  adverse  pressure  rise 
(due  to  p„2  increasing  while  Pc\  decreases)  and  decrease  in 
streamwise  momentum  flux  (due  to  decreasing)  with  in¬ 
creasing  Mach  number. 

The  ability  of  the  model  to  predict  the  onset  of  shock- 
induced  breakdown  is  further  evaluated  in  Fig.  3,  where  data 
from  Metwally  era/."  are  plotted  (obtained  from  Rg.  6  of 
their  paper).  The  “strong  interactions’’  observed  experimen¬ 
tally  ate  seen  to  lie  in  the  region  where  the  model  predicts 
breakdown,  while  the  “weak  interaction’’  regions  lie  in  flie 
predicted  region  of  non-breakdown.  Note  that  the  curve  of 
Ten,  in  Fig.  3  assumes  uniform  axial  velocity.  Metwally 
et  a/."  point  out  that  the  Mach  3  and  Mach  33  vortices  had 
noticeable  deficit  in  centerline  velocity  for  the  breakdown 
cases.  As  will  be  seen  in  Section  DIB,  the  critical  swirl 
number  is  predicted  to  decrease  as  the  centerline  velocity 
decreases;  i.e.,  the  filled  symbols  for  the  Mach  3  and  Mach 
3.5  cases  would  move  further  into  the  breakdown  region  if 
the  deficit  in  centerline  velocity  were  accounted  for  in  Hg.  3. 


B.  Non-uniform  axial  velocity 

The  influence  of  an  excess/deficit  in  the  centerline  axial 
velocity  on  the  critical  swirl  number  is  next  considered.  For 


1 
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FIG.  2.  ^mparison  of  predicted  critical  swirl  number  to  experiment  and  FIG.  3.  Evaluation  of  model  in  predicting  the  onset  of  shock-induced  vortex 

computotton  of  shock-indu(*d  vortex  breakdown.  —  (Prediction:  unifonn  breakdown.  —  (Predicted  [Eq.  (17)]:  uniform  stagnation  temperature), 

stagnation  temperature).  —  (prediction:  uniform  entropy),  •  •  (experiment— Ref.  11:  breakdown),  O  (experiment— Itef.  11:  no  break- 

(compuution-Ref.  17),  X  (experiment— Ref.  10).  down).  • 
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FIG,  4.  Influence  of  axial  velocity  on  the  onset  of  vortex  breakdown  in¬ 
duced  by  a  shock.  —  (Af//C/»,  =  -0.5).  —  (AC//{/*|=  “0.25),  — 
(At//(/ac|=0),  —  (Af//f/«,=0.5),  (Al//I/»|  =  1). 


convenience,  results  are  shown  only  for  the  uniform  entropy 
vortex.  The  breakdown  criterion  [Eq.  (12)]  may  be  divided 
through  by  /7x|  to  yield  the  following  non-dimensional  cri¬ 
terion: 

'I  AC/\^ 

Po=2-Pci  =  7PciMli^l  +  -jj^J  ■  (19) 

Substituting  for  and  from  Eqs.  7,  we  get  the 
following  equation  for  Fcni  as  a  function  of  Afooi 
Af//t/aoi: 

=yA#2,[  1+ ^j\i  (20) 

The  Newton-Raphson  method  was  used  to  solve  the 
above  equation  for  Tent,  after  expressing  p„2  in  terms  of 
Ma,\.  Figure  4  shows  the  variation  of  the  critical  swirl  num¬ 
ber  with  Mach  number  for  different  values  of  A  f//t/oo  i .  Note 
that  Atf>0  corresponds  to  a  jet-like  axial  velocity  profile  of 
the  upstream  vortex  while  Af/<0  corresponds  to  a  wake¬ 
like  profile.  The  predicted  results  show  a  strong  sensitivity  to 
the  excess/deficit  in  centerline  axial  velocity.  Jet-like  profiles 
of  the  axial  velocity  are  observed  to  delay  breakdown,  while 
a  wake-like  profile  makes  the  vortex  more  susceptible  to 
breakdown.  The  same  trend  is  known  to  apply  in  the  break- 
doivn  of  an  incompressible  vortex,  where  axial  blowing  is 
often  used  to  alleviate  the  breakdown.^ 

Figure  4  shows  that  for  vortices  with  a  wake-like  axial 
velocity,  the  critical  swirl  number  becomes  zero  at  a  finite 
Mach  number;  i.e.,  breakdown  is  predicted  at  and  beyond 
this  cut-off  Mach  number,  even  in  the  absence  of  swirl .  This 
result  may  be  explained  as  follows.  In  the  absence  of  swirl, 
the  “vortex”  reduces  to  an  axisytmnetric  wake  (or  jet).  This 
wake  (or  jet)  can  undergo  reverse  flow  accompanied  by  ra¬ 
dial  outflow  upon  experiencing  a  strong  enough  adverse 
pressure  gradient.  We  have  assumed  that  breakdown  occurs 
when  the  adverse  pressure  rise  at  the  vortex  centerline  ex- 
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ceeds  the  centerline  axial  momentum  flux.  We  noted  that  on 
account  of  its  rotation,  the  centerline  pressure  rise, 
Foo2~Pci  is  greater  than  the  free-stream  rise,  p„2-pm\- 
Also,  rotation  results  in  the  centerline  density  (Pd)  being 
lower  than  the  ftee-stream  density.  As  a  result,  the  centerline 
momentum  flux,  Pc\U\x  is  less  than  its  value  computed  us¬ 
ing  the  free-stream  density.  Thus,  swirl  “amplifies”  (in  the 
terminology  of  Hall')  the  adverse  pressure  rise  experienced 
by  the  vortex  while  suppressing  the  axial  momentum  flux. 
Both  factors  make  the  vortex  more  susceptible  to  breakdown. 
This  implies  that  if  the  free-stream  pressure  rise  exceeds  the 
axial  momentum  flux  computed  using  the  free-stream  den¬ 
sity,  then  the  presence  of  swirl  is  not  needed  for  “break¬ 
down.”  The  flow  at  and  above  the  cut-off  Mach  number 
corresponds  to  this  scenario.  The  cut-off  Mach  number  (de¬ 
noted  by  Afeut)  can  therefore  be  predicted  by  the  following 
criterion: 

Poo2~Pool  “P*l  (21) 

which  yields, 

,  /  Af/\2 

Px2-l  =  rMHl  +  j^j  .  (22) 

Substituting  for  p„2  from  Eq.  (15),  we  get  the  following 

equation  for  the  cut-off  Mach  number  in  terms  of  the  veloc¬ 
ity  excess/deficit: 

2y  ,  2  I 

+  (23) 

which  may  be  rearranged  to  obtain  the  following  expression 
for  the  cut-off  Mach  number: 


C.  Breakdown  induced  by  an  oblique  shock  wave 

If  the  shock  wave  were  oblique,  the  onset  of  breakdown 
would  be  expected  to  depend  on  the  oblique  angle.  Although 
the  interaction  of  an  oblique  shock  with  an  axisymmetric 
vortex  is  not  axisymmetric,  it  is  envisioned  that  the  onset  of 
breakdown  can  be  predicted  by  extending  the  arguments  of 
the  previous  section.  Reiterating  the  criterion  for  breakdown 
for  uniform  axial  velocity,  we  require  that 
^oo2‘“^ci“  Th^  influence  of  shock  obliquity  is 

modeled  as  follows.  The  properties  of  the  upstream  vortex 
ip  cl  y  Pci)  depend  solely  upon  the  free-stream  Mach  number 
and  swirl  number.  However  the  pressure  behind  the  shock 
(^002)  is  determined  by  the  normal  Mach  number,  Afocisina 
(a  denotes  the  angle  the  shock  makes  with  the  streamwise 
direction).  Replacing  Afeej  in  Eq.  (15)  by  Af^isina  to  obtain 
p^2  ^d  substituting  as  before  for  Pc\  and  p^y  yields  the 
following  expressions  for  the  critical  swirl  number. 
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Uniform  stagnation  temperature  vortex: 


1  /  2  (  /  i 


sin 


VFTWy 


Uniform  entropy  vortex: 

‘  +  sina]2-l)-[l-(y-l)r^^2jr/(y-i)=yM^,[l_(y-l)r2^^2,]i/(r-i). 


(25) 

(26) 


It  is  readily  seen  that  for  the  same  upstream  Mach  num- 
her,  r  crit  is  predicted  to  increase  as  the  shock  becomes  in¬ 
creasingly  oblique.  This  prediction  may  be  explained  by  not¬ 
ing  that  the  pressure  rise  across  an  oblique  shock  is  lower 
than  that  for  a  normal  shock  at  the  same  Mach  number.  As  a 
result,  the  adverse  pressure  rise  that  the  vortex  experiences  is 
smaller,  thereby  delaying  the  onset  of  breakdown. 


IV.  SHOCK-FREE  BREAKDOWN  OF  A 
COMPRESSIBLE  VORTEX 

Section  HI  discussed  vortex  breakdown  induced  by  a 
shock  wave.  The  breakdown  of  a  free  axisymmetric  vortex, 
i.e.  breakdown  in  the  absence  of  an  externally  imposed  pres¬ 
sure  gradient,  is  considered  in  this  section.  Incompressible 
streamwise  vortices  at  sufficiently  high  swirl  number  are 
known  to  break  down,  even  in  the  absence  of  an  externally 
applied  adverse  pressure  gradient.  It  is  to  be  expected  that 
their  high-speed  counterparts  would  exhibit  similar  behavior. 
The  critical  swirl  number  in  high-speed  flow  would  be  a 
function  of  the  Mach  number.  This  section  derives  an  ex¬ 
pression  for  the  critical  swirl  number  in  terms  of  the  firee- 
stream  Mach  number;  i.e.,  we  consider  the  influence  of  com¬ 
pressibility  on  the  breakdown  of  a  free  vortex.  The 
arguments  used  are  identical  to  those  in  breakdown  induced 
by  a  shock.  The  only  difference  is  that  while  the  adverse 
pressure  rise  was  set  equal  to  Poo2~Pci  for  shock-induced 
breakdown,  it  is  set  equal  to  Poo|-“Pci  for  the  shock-free 
breakdown.  The  rationale  for  this  assumption  is  that  in  the 
absence  of  the  shock,  the  vortex  discharges  into  the  atmo¬ 
sphere.  As  a  result,  the  vortex  sees  a  pressure  equal  to  p^x 
ahead  of  it,  as  well  as  in  the  free-stream.  The  difference 
between  atmospheric  pressure  (poci),  and  the  pressure  at  the 
vortex  centerline  (p^j)  provides  the  adverse  pressure  rise 
that  causes  breakdown.  Breakdown  of  the  vortex  is  therefore 
assumed  to  occur  when 

P^\-Pc\^Pc\U]x.  (27) 

The  criterion  for  shock-free  breakdown  is  therefore  given  by, 

l“Pci  =  yPciA^xi  (28) 

which  is  identical  to  the  expression  obtained  when  ^„2  is  set 
to  1  in  Eq.  (13).  The  corresponding  expressions  for  the  criti¬ 
cal  swirl  number  are  given  below. 


Uniform  stagnation  temperature  vortex: 

1  /  2  [  1  i 


(29) 


Uniform  entropy  vortex: 

=yM^,[i-(y-i)rc^ri,wii]‘'‘’'~‘>.  (30) 

Figure  5  shows  the  predicted  values  of  the  critical  swirl 
number  as  a  function  of  the  free-stream  Mach  number.  Also 
shown  (for  supersonic  flow)  are  the  values  obtained  for 
breakdown  induced  by  a  shock  wave  at  the  same  Mach  num¬ 
ber.  Compressibility  is  seen  to  make  the  vortex  more  suscep¬ 
tible  to  breakdown.  A  similar  trend  was  noted  by  Keller.^^ 
This  trend  may  be  explained  by  noting  [Eqs.  (7)  and  (10)] 
that  increase  in  the  free-stream  Mach  number  decreases  the 
centerline  pressure  and  density,  thereby  increasing  the  ad¬ 
verse  pressure  rise  while  decreasing  the  axial  momentum 
flux.  Tlie  predicted  values  of  F  ^nt  t*t  the  absence  of  the  shock 
are  seen  to  be  greater  than  those  predicted  for  shock-induced 
breakdown.  This  trend  can  be  explained  by  noting  that  the 
pressure  rise  across  the  shock  wave  produces  a  larger  ad¬ 
verse  pressure  rise  for  the  same  upstream  momentum  flux. 


FIG.  5.  Predicted  critical  swirl  number  for  shock-free  vortex  breakdown 
compared  to  the  prediction  for  shock-induced  breakdown.  —  (Shock-free: 
uniform  stagnation  temperature),  —  (shock-free:  uniform  entropy), 
(shock-induced:  uniform  stagnation  temperature),  — -  (shock-induced:  uni¬ 
form  entropy). 
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TABLE  L  Prediction  of  critical  swirl  number  for  incompressible  vortex 
breakdown  compared  to  other  approaches.  All  data  other  than  the  present 
lepnxhiced  from  review  article  by  Dclery  (Ref.  6). 


1.0 


Quasi-cylindrical 

1.41 

Axisynunetric  N*S 

1.35 

Bossel 

1.12 

Squire 

1.4 

Benjamin 

1.4 

Norn,  simulation 

1.28 

Spall  et  al. 

1.37 

Present 

1.4 

V.  INCOMPRESSIBLE  VORTEX  BREAKDOWN 

Hgure  5  shows  that  as  tends  towards  0,  Tent  tends 
towards  1.  An  incompressible  vortex  in  the  absence  of  exter¬ 
nally  imposed  adverse  pressure  gradients,  is  therefore  pre¬ 
dicted  to  undergo  breakdown  at  a  critical  swirl  number  of 
one.  The  same  result  can  of  course  be  derived,  by  setting 
p—poo  in  the  radial  momentum  equation  and  integrating  to 
obtain  the  centeriine  pressure  iPcx—Pooi-^ PoaV%n)*  which  is 
then  substituted  into  the  breakdown  criterion  [^.  (27)].  In  a 
recent  review  article,  Delery^  documents  (Section  3.4.5  of 
his  paper)  critical  swirl  numbers  for  incompressible  vortex 
breakdown  as  predicted  by  different  theories.  He  considers  a 
Burgers  vortex,  and  defines  a  swirl  parameter  S  as 


where  the  variables  C  and  denote  die  circulation  and  core 
radius  respectively.  For  a  Burgers  vortex,  the  swirl  velocity 
is  given  by  (Eq.  1  in  Delery's^  paper) 


1.256(r/r^)2j 


(32) 


This  implies  that  the  swirl  parameter  S  is  related  to  the  swirl 
number  F  by. 


L398r. 


(33) 


Thus  rcrit==l  corresponds  to  5crit=1.398«=»1.4.  We  repro¬ 
duce  in  Table  I,  from  Delery’s^  paper,  the  critical  swirl  num¬ 
bers  predicted  by  different  sqiproaches.  Most  approaches  are 
seen  to  predict  values  very  close  to  that  predicted  by  our 
simple  criterion. 


VI.  A  “UNIVERSAL”  BREAKDOWN  MAP 

The  preceding  sections  presented  results  for  the  onset  of 
vortex  breakdown  by  plotting  the  critical  swirl  number  as  a 
function  of  Mach  number.  The  curve  Tcnt—^cn^Mooi)  de¬ 
fined  the  boundary  between  the  regimes  of  breakdown  and 
non-breakdown.  However,  it  is  clear  that  the  critical  swirl 
number  is  not  universal  (as  also  observed  by  Delery^).  For 
example.  Section  m  B  (Fig.  4)  showed  that  F crit  depended  on 
the  velocity  excess/deficit  at  the  centerline.  If  the  breakdown 
were  precipitated  by  an  oblique  shock  wave  as  opposed  to  a 
normal  shock,  then  F^it  was  noted  to  depend  on  the  inclina¬ 


FIG.  6.  Evaluation  of  the  proposed  breakdown  map  in  predicting  tiie  onset 
of  vortex  breakdown.  #  (Experiment:  breakdown),  O  (experiment:  no 
breakdown). 


tion  angle  of  the  shock.  Similarly,  if  the  breakdown  were  that 
of  a  &ee  vortex  instead  of  being  shock  induced,  yet  another 
curve  for  the  critical  swirl  number  was  obtained. 

In  this  section,  we  propose  a  breakdown  that  allows 
a  common  breakdown  boundary  to  be  defined  for  all  of  the 
above  mentioned  problems.  The  proposed  map  is  based  on 
the  breakdown  criterion  that  was  proposed  in  Section  n  B; 
i.e., 

P'»2-Pci^PciUli.  (34) 

Recall  that  die  same  criterion  with  t^ipropriately  defined, 
was  iqiplied  to  all  the  breakdown  problems  discussed  in  this 
paper.  This  suggests  that  a  plot  of  Paa2~Pci  against  PciUh 
could  be  used  to  map  the  onset  of  vortex  breakdown.  The 
proposed  map  could  even  be  used  for  incompressible  vortex 
breakdown,  and  would  be  expected  to  adequately  represent 
the  onset  of  breakdown  induced  by  pressure  gradients  acting 
over  distances  that  are  small  as  compared  to  a  characteristic 
length  scale  of  the  vortex.  The  curve  Poa2~Pci~Pcif^ci 
45°  line)  would  act  as  die  boundary  between  the  breakdown 
and  non-breakdown  regimes.  Note  that  the  proposed  map 
does  not  require  any  additional  data  to  be  measur^.  Experi¬ 
mental  informadon  cm  parameters  such  as  T,^UIU^,M^ 
and  shock  angle  could  be  used  to  obtain  both  the  pressure 
rise  and  the  axial  momentum  flux  using  the  equations  in 
Section  n  A.  The  proposed  map  is  illustrated  in  Fig.  6.  Note 
that  the  pressure  rise  and  momentum  flux  are  non- 
dimensionalized  by  Pa|C/«|  to  allow  incompressible  data  to 
be  plotted.  Data  firom  Metwally  (the  same  data  shown  in 
Fig.  3)  are  also  shown.  The  data  finom  Hg.  3  are  combined 
with  Eq.  10  to  determine  the  pressure  rise  and  axial  momen¬ 
tum  flux.  The  breakdown  and  non-breakdown  cases  ate  seen 
to  be  appropriately  delineated. 

VII.  SUMMARY 

A  simple  inviscid  model  was  proposed  to  predict  flie 
onset  of  breakdown  in  an  axisynunetric  vortex.  Three  prob¬ 
lems  were  considered:  the  shock-induced  breakdown  of  a 
compressible  vortex,  the  breakdown  of  a  fiee  compressible 
vortex,  and  the  breakdown  of  a  free  incompressible  vortex. 
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The  same  physical  reasoning  was  used  to  predict  the  onset  of 
breakdown  in  all  three  problems.  It  was  hypothesized  that 
breakdown  is  the  result  of  the  competing  effects  of  adverse 
pressure  rise  and  streamwise  momentum  flux  at  the  vortex 
centerline.  Breakdown  was  assumed  to  occur  if  the  pressure 
rise  exceeded  the  axial  momentum  flux.  A  formula  with  no 
adjustable  constants  was  derived  for  the  critical  swirl  number 
in  all  three  problems.  The  dependence  of  the  critical  swirl 
number  on  parameters  such  as  upstream  Mach  number, 
excess/deficit  in  centerline  axial  velocity,  and  shock  oblique 
angle  was  explored.  The  predictions  for  the  onset  of  shock- 
induced  breakdown  and  free  incompressible  breakdown  were 
compared  to  experiment  and  computation,  and  good  agree¬ 
ment  was  observed.  Finally,  a  new  breakdown  map  was  pro¬ 
posed  as  an  alternative  to  the  map  of  critical  swirl  number 
against  free-stream  Mach  number.  The  new  map  was  based 
on  the  observation  that  the  same  breakdown  criterion  was 
used  in  all  the  problems  considered  in  this  paper.  To  deter¬ 
mine  the  onset  of  breakdown,  it  was  suggested  that  the  ad¬ 
verse  pressure  rise  at  the  vortex  centerline,  be  plotted  against 
the  axial  momentum  flux.  The  proposed  map  allows  the  si¬ 
multaneous  comparison  of  data  from  flows  ranging  from  in¬ 
compressible  breakdown  to  breakdown  induced  by  a  shock 
wave. 
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ABSTRACT 

This  paper  presents  a  fenoily  of  finite  difierence  schemes  for  the  first  and  second 
derivatives  of  smooth  functions.  The  schemes  are  Hermitian  and  symmetric,  and  may  be 
considered  a  more  general  version  of  the  standard  compact  (Fade)  schemes  discussed  by 
Lele  [1].  They  are  different  fi:om  the  standard  Fade  schemes,  in  that  the  first  and  sec¬ 
ond  derivatives  are  evaluated  simultemeously.  For  the  same  stencil  width,  the  proposed 
schemes  are  two  orders  higher  in  accuracy,  and  have  significantly  better  spectral  represen¬ 
tation.  Eigenvalue  analysis,  and  nmnerical  solutions  of  the  one-dimensional  wave  equation 
are  used  to  demonstrate  the  numericed  stability  of  the  schemes.  The  computational  cost 
of  computing  both  derivatives  is  assessed,  and  shown  to  be  essentially  the  same  as  the 
standard  Fade  schemes.  The  proposed  schemes  appear  to  be  attractive  alternatives  to  the 
standard  Fade  schemes  for  computations  of  the  Navier  Stokes  equations. 

1.  INTRODUCTION 

Fluid  fiows  in  the  transitional  and  ttirbulent  regimes  possess  a  wide  range  of  length 
and  time  scales.  The  numerical  computation  of  these  flows  therefore  requires  numerical 
methods  that  can  accurately  represent  the  entire,  or  at  least  a  significant  portion,  of  this 
range  of  scales.  The  length  scales  that  are  resolved  by  a  computation  are  determined  by  the 
resolution;  the  acctiracy  with  which  these  scales  are  represented  depends  upon  the  numer¬ 
ical  scheme.  Fourier  anedysis  (see  e.g.  [2])  describes  both,  the  range  of  scales  present,  and 
the  accuracy  with  which  they  are  computed  (exactly  for  problems  with  periodic  boundary 
conditions,  and  in  a  WKB  sense  for  more  general  problems).  Such  analysis  of  finite  differ¬ 
ence  schemes  (see  e.g.  Fig.  1  in  [1])  shows  that  the  error  in  computing  the  first  and  second 
derivatives  can  be  quite  large  for  the  smaller  scales.  This  small  scale  inaccuracy  becomes 
increasingly  important  as  the  energy  in  the  small  scales  becomes  increasingly  comparable 
to  that  of  the  large  scales;  t.e.,  as  the  spectrum  becomes  increasingly  ‘flat’.  This  situ¬ 
ation  is  coinmordy  encountered  in  computations,  particrdarly  large-eddy  simulations,  of 
high  Reynolds  number  turbulence.  As  shown  by  Kravchenko  and  Moin  [3]  the  inaccmate 
numerical  representation  of  the  small  scales  in  these  large-eddy  simulations  can  result  in 
the  numerical  error  overwhelming  the  contribution  of  the  subgrid-scale  model. 

Finite  difference  schemes  may  be  classified  as  ‘explicit’  or  ‘implicit’.  Explicit  schemes 
express  the  nodal  derivatives  as  an  explicit  weighted  sum  of  the  nodal  values  of  the  function, 
e.ff.,  f'i  =  ifi+i  -  /i-i)/2h,  and  //'  =  {fi+i  -  2/,-  +  fi-i)lh'^.  Throughout  this  paper,  f, 
and  fi  denote  the  values  of  the  function  and  its  derivative  respectively,  at  the  node 
=  Xi,  and  h  denotes  the  imiform  mesh  spacing.  By  comparison,  implicit  (compact) 


X 
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equate  a  weighted  sum  of  the  nodal  derivatives  to  a  weighted  sum  of  the  function,  e.g., 
fU  +^fi  +  fUi  =  3(/.+i  -fi-i)/h,  and  + 10//'  +  /^;,  =  12(/i+i  -2fi+ fi_i)/h^.  It 

is  well  known  [1,4,5]  that  implicit  schemes  have  better  small  scale  accuracy,  than  explicit 
schemes  with  the  same  stencil  width.  This  increase  in  accuracy  is  achieved  at  the  cost 
of  inverting  a  banded  (usually  tridiagonal)  matrix  to  obtain  the  nodal  derivatives.  Since 
tridiagonal  matrices  can  be  inverted  quite  efficiently  [6],  the  implicit  schemes  are  extremely 
attractive  when  explicit  time  advancement  schemes  are  used.  The  most  popular  of  the 
implicit  schemes  (also  called  Fade  schemes  due  to  their  derivation  from  Fade  approximants) 
appear  to  be  the  symmetric  fotirth  and  sixth  order  versions  (see  e.g.  [1]).  There  have  been 
several  recent  computations  of  transitional  boundary  layers  [7-10],  turbulent  flows  [11-13] 
and  flow-generated  noise  [14-15]  that  have  used  the  Fade  schemes  to  evaluate  the  spatial 
derivatives.  The  standard  Fade  schemes  are  symmetric  and  therefore  non-dissipative;  a 
non-symmetric  compact  scheme  was  recently  developed  by  Adams  and  Shariff'  [16]. 

This  paper  presents  a  related  family  of  finite  difierence  schemes  for  the  spatial  deriva¬ 
tives  in  the  Navier  Stokes  equations.  The  proposed  schemes  are  more  accurate  than  the 
standard  Fade  schemes,  while  incurring  essentially  the  same  computational  cost.  They  are 
based  on  Hermite  interpolation,  and  may  be  considered  a  more  general  version  of  the  stan¬ 
dard  Fade  schemes  described  in  [1].  For  the  same  stencil  width  as  the  Fade  schemes,  the 
proposed  schemes  have  higher  order  of  accuracy  and  better  spectral  representation.  This 
is  achieved  by  simultaneously  solving  for  the  first  and  second  derivatives.  When  defined 
on  a  uniform  mesht  ,  the  schemes  are  of  the  form. 


Ol/i-l-|-ao//+a2/i+i+h(&l  /<1i+&o/"+62/h.i)  =  ^(ci/,_2+C2/i_l-|-Co/i-l-C3/,+iH-C4/,+2) 

.r  . 

Note  that  the  above  expression  differs  from  the  standard  Fade  schemes,  in  that  the 
left-hand  side  contains  a  linear  combination  of  the  first  and  second  derivatives.  The  stencil 
and  the  coefficients  are  restricted  to  be  symmetric  in  this  paper.  The  resulting  schemes 
are  therefore  non-dissipative.  The  width  of  the  stencil  is  taken  to  be  three  on  the  left-hand 
side  and  five  on  the  right.  This  corresponds  to  the  stencil  width  of  the  popular  sixth-order 
Fade  scheme. 

The  motivation  to  formulate  schemes  that  simultaneously  evaluate  both  derivatives  is 
provided  by  the  Navier  Stokes  equations  requiring  both  derivatives  of  most  variables.  Con¬ 
sider  for  example  the  one-dimensional  compressible  equations  in  primitive  form  (extension 
to  multiple  dimensions  is  straightforward).  We  have: 


dp  dp  du 
dt  ^  dx  ^  dx’ 


(2a) 


t  This  paper  develops  the  schemes  on  uniform  meshes.  It  is  assumed  that  computations  with 
non-uniform  grids  can  define  analytical  mappings  between  the  non-uniform  grid  and  a  corresponding 
uniform  grid.  The  metrics  of  the  mapping  may  then  be  used  to  relate  the  derivatives  on  the  uniform 
grid  to  those  on  the  non-uniform  grid. 
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{dn  du\  ^dp  „ 

^dx  ^dx  ZdxdTdx, 


0-(iaSiL\ 

Ox  V  3  ^  Ox  } 

r&T  &T\  du  4  fduV  .a^T  dh  ( dTV 

"^dT  (ax)  ;  ^ 

The  variables  p,u  and  T  denote  the  density,  velocity  and  temperature  respectively,  while 
R,p,k  and  C„  denote  the  specific  gas  constant,  dynamic  viscosity,  thermal  conductivity 
and  specific  heat  at  constant  volume.  Note  that  the  viscous  terms  axe  expanded  prior  to 
their  evaluation.  This  is  because  direct  evaluation  of  the  second  derivatives  is  sigmficantly 
more  accurate  than  two  applications  of  a  first  derivative  operator.  Gquation  2  shows  that 
the  following  spatial  derivatives  need  to  be  evaluated: 

^  ^  ^  ^  ^ 

dx  ’  dx"^ '  dx  ’  dx^  ’  dx 

Thus,  a  scheme  that  simultaneously  evaluates  both  derivatives  would  only  be  performing 
one  unnecessary  evaluation 

Next,  consider  the  conservative  form  of  the  equations.  The  viscous  terms  are  still 
evaluated  in  their  non-conservative  form,  for  the  reasons  given  above.  We  have: 

^  +  ^(P«)=0-  (3«) 

9  ,  _  9  ,  .  9p_A__di^u  4^^^  .... 


dx  dx^  ZdxdTdx 


dEt  ^ 
dt  ^  dx 


d  /'  /4  d^u  4  5m  d/i  5r\  4  { 9u\ 

{E,n)  +  ^  (pm)  =m  -h  +  3^ 

d^T 

■^^5x2  dT  V5xy  ‘ 


Equation  3  requires  the  following  spatial  derivatives  to  be  obtained: 

a,  ,  a  ,  ’  ,  9p  9u  a^u  ar  a=r  9 


,  „  dp  du  a^u  ar  a=r  d  .  ,  9,  , 

9i’  S’  s?’  S’  s?’  s‘^‘“>’  s(p“) 


As  one  might  expect,  the  conservative  formulation  requires  fewer  simultaneous  derivative 
evaluations.  However,  if  the  chain  rule  is  invoked  as  follows,  then  a  formulation  that 
evaluates  both  derivatives  is  still  attractive.  First  evaluate  (simultaneously) 


52  .  d  f  2\ 

ipu). 


.  o\  9p  d^p  9  .  52 

£(pm),  and  ^(pM). 
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The  chain  rule  may  then  be  used  to  obtain  du/dx^  d^u/dx^,  dpjdx  and  d^p/dx^.  The 
equation  of  state  and  the  chain  rule  then  yield  dTfdx  and  df^T/dx^.  In  this  manner,  a  total 
of  only  ten  derivative  evaluations  are  performed  for  the  nine  derivatives  that  are  needed. 
The  increase  in  accuracy  that  is  obtained  by  the  simultaneous  evaluation  of  derivatives 
will  be  seen  to  maJce  this  additional  derivative  evaluation  worthwhile. 

For  the  same  stencil  width,  the  standeird  Fade  schemes  are  two  orders  higher  in  accu¬ 
racy  amd  have  better  spectral  representation  than  the  corresponding  symmetric,  explicit 
schemes.  The  implicit  relation  between  the  derivatives  in  the  Fade  schemes  yields  addi¬ 
tional  degrees  of  freedom  that  allow  higher  accuracy  to  be  achieved.  It  is  therefore  to 
be  expected  that  including  the  second  derivatives  in  the  implicit  expression  would  further 
increase  the  degrees  of  freedom,  and  thereby  the  accuracy  that  can  be  obtained.  Hermitian 
expressions  involving  functions  and  their  first,  and  higher  derivatives  have  been  suggested 
in  the  literature  {e.g.  [4],  sections  2.4,  2.5).  Feyret  and  Taylor  ([17],  section  2.5.1)  and 
Hirsch  ([18],  section  4.3)  discuss  a  symmetric  version  of  equation  1  on  a  three  point  stencil. 
However,  the  development  was  not  completed  to  a  point  where  the  resulting  schemes  could 
be  used  for  solving  partial  differential  equations. 

The  objective  of  this  paper  is  to  develop  this  family  of  schemes,  and  assess  their 
potential  for  computations  of  the  Navier  Stokes  equations.  The  schemes  will  be  referred  to 
as  the  ‘coupled-derivative’,  or  ‘C-D’  schemes  to  distinguish  them  from  the  standard  Fade 
sdiemes.  The  paper  is  organized  as  follows.  Section  2  describes  the  interior  schemes  that 
may  be  obtained  from  equation  1.  Fourier  analysis  is  then  used  in  section  3  to  perform 
a  detailed  comparison  between  the  proposed  schemes  and  the  standard  Fade  schemes. 
The  restrictions  imposed  by  numerical  (Cauchy)  stability  are  then  discussed  in  section  4. 
Section  5  presents  appropriate  boundary  closures  for  the  interior  scheme,  and  evaluates 
the  stability  of  the  complete  scheme.  The  computational  cost  of  the  proposed  schemes  is 
evaluated  in  section  6,  and  compared  to  that  of  the  standard  Fade  schemes.  The  paper  is 
concluded  with  a  brief  summary  in  section  7. 

2.  THE  INTERIOR  SCHEME 

The  interior  scheme  is  of  the  form  given  by  equation  1.  Simultaneous  solving  for  f- 
and  /",  implies  that  the  n\imber  of  unknowns  is  equal  to  2N.  A  total  of  2N  equations 
are  therefore  needed  to  close  the  system.  Equation  1  may  be  used  to  derive  two  linearly 
independent  equations  at  each  node.  This  is  done  as  follows.  Both  sides  of  equation  1 
are  first  expanded  in  a  Taylor  series.  The  resulting  coefficients  are  then  matched,  such 
that  equation  1  maintains  a  certain  order  of  accuracy.  Note  that  equation  1  has  eleven 
coeflS-cients,  of  which  one  is  eurbitrary,  i.  e.,  equation  1  may  be  divided  through  by  one  of 
the  constants,  without  loss  of  generality.  A  convenient  choice  of  the  normalizing  constant, 
is  either  of  oq  or  feo-  It  will  be  seen  that  the  equation  obtained  by  setting  oq  equal  to  1,  is 
linearly  independent  of  the  equation  obtained  when  bo  is  set  equal  to  1.  The  two  equations 
may  therefore  be  applied  at  each  node,  and  the  resulting  system  of  2N  equations  solved 
for  the  nodal  values  of  the  first  and  second  derivative.  The  process  of  obtaining  the  two 
equations  is  outlined  in  sections  2.1  and  2.2. 
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fr 

0 

0 
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2h»  (01/8! -b  62/7!) 

2/1®  (2®  C4  -b  C3)  /9! 

TABLE  1:  Taylor  table  for  Oq  =  1. 


2.1:  First  equation  (ao  =  \) 

Consider  first  the  case  where  ao  =  1.  The  symmetry  of  the  schemes  requires  that 
oi  =  02,  6i  =  —62,  Cl  =  — C4  and  C2  =  —c$.  Equation  1  therefore  reduces  to  the  form: 

Co/,+C3(/,--(.i— /f_i)+C4(/,+2— /i-2)j 

(4 

Expanding  both  sides  of  equation  4  in  a  Taylor  series  and  collecting  terms  of  the  same 
order  yields  Table  1.  Note  that  ‘LHS’  and  ‘RHS’  denote  the  coeflScients  of  /*  on  the  left 
and  right-hand  sides  respectively  of  equation  4. 

The  Taylor  table  shows  that  60  =  co  =  0.  This  leaves  four  undetermined  constemts 
(oi ,  62 ,  C3  and  C4).  Expressions  for  these  constants  may  be  obtained  by  matching  the  terms 
in  the  Taylor  table.  Schemes  of  order  ranging  from  two  through  eight  may  be  obtained  by 
solving  the  resulting  set  of  equations.  The  coefficients  and  the  resulting  orders  are  listed 
below. 

Second  order 

Matching  terms  up  to  /"  yields, 


01  =  — -  -I-  C3  +  2c4,  &2  arbitrary. 

St 

The  resulting  leading  order  error  is  equal  to  (3  —  1262  —  4c3  +  4Ci)hr  /6. 


(5u) 
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Fourth  order 

Matching  terms  up  to  /?"  yields, 

«i  =  -^  +  C3  +  2c4,  62  =  ^[3  -  4(c3  -  C4)].  (56) 

The  resulting  leading  order  error  is  given  by(— 15  +  I6C3  +  92c4)/i^/"/360.  Note  that 
C4  =  0,  C3  =  3/4  yields  the  standard  fourth  order  Fade  scheme  for  the  first  derivative. 

Sixth  order 

Matching  terms  up  to  yields, 


7  15  .  1  r  ,  N  15  23 

ai  =  —  -  — C4,  62  =  jg(-l  +  36c4),  =  Yg  -  ^C4. 


(5c) 


The  resulting  leading  order  error  is  equal  to  (l/5040+3c4/140)6®/?’”.  Note  that  C4  =  1/36 
yields  the  standard  sixth  order  Fade  scheme  for  the  first  derivative. 

Eighth  order 

Matching  terms  up  to  /"”*  yields. 
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107 

108’ 


C4 


1 

108  ■ 


The  error  to  leading  order  is  equal  to  — h®/?^/90720. 

Table  1  shows  that  60  is  equal  to  zero  when  Oo  is  set  equal  to  one.  The  above  expres¬ 
sions  may  therefore  be  considered  expressions  for  the  nodal  values  of  the  first  derivative. 
It  also  implies  that  if  instead  of  setting  oq  equal  to  one,  we  set  60  equal  to  one,  we  would 
obtain  an  equation  that  would  be  linearly  independent.  The  equation  thus  derived  could 
be  considered  an  expression  for  the  second  derivative.  This  equation  is  obtained  below. 


2.2:  Second  equation  (bo  =  \) 

Consider  the  case  where  60  =  1.  Note  that  a  tilde  is  used  above  the  constants  to 
indicate  their  difference  frorn  the jionst ants  obtained  when  oq  =  1;  e.g.,  h\  is  replaced  by 
61.  Symmetry  requires  that  61  =62,  Ci  =  04,  c%  =  C3  and  di  =  —02.  Equation  1  therefore 
becomes: 

Cl  (/»-2 + /*+2 ) + C2  (/i_l -f/i+i ) -bco/i 

Expanding  both  sides  of  the  above  equation  in  a  Taylor  series  and  collecting  terms  of  the 
same  order  yields  the  Taylor  table  2. 

Table  2  shows  that  cq  is  required  to  be  zero  if  60  is  equal  to  one.  The  resulting  equation 
may  therefore  be  considered^an  expression  for  the  second  derivative.  We  have  five  unknown 
constants  (cq  ,  ci ,  C2 , 02  and  64 ).  These  constants  may  be  obtained  by  matching  the  terms  in 
the  above  Taylor  table,  and  solving  the  resulting  equations.  Expressions  of  varying  order 


aof'+MfUi-fU)+h{b^  //L4+/r+64/r;j  =  ± 
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®  TABLE  2:  Taylor  table  obtained  for  6o  =  1* 


axe  obtained,  depending  upon  the  number  of  equations  matched.  At  first  glance,  it  appears 
that  the  order  of  accuracy  obtained,  ranges  from  three  through  nine.  By  comparison,  the 
expressions  obtained  when  Oq  was  equal  to  1  ranged  from  second  through  eighth  order. 
However,  note  that  the  nodal  second  derivatives  in  equation  1  are  premultiplied  by  h. 
Equation  1  (and  therefore  the  terms  in  the  Taylor  table)  needs  to  be  divided  through  by  h, 
to  consider  it  an  expression  for  the  second  derivatives.  This  process  will  yield  expressions 
for  the  second  derivative,  ranging  in  order  from  two  through  eight.  The  values  for  the 
constants  and  the  corresponding  orders  are  given  below. 


Second  order 


Matching  terms  up  to  /"  yields. 


Co  =  -2(ci  +  ca),  0.2  =  2 ^2)- 

The  resulting  leading  order  error  is  (2  —  861  +  8ci  —  C2)h^/,*’'/12. 

Fourth  order 

Matching  terms  up  to  yields, 

3  5..  r  1  ~  C2 

Co  = -2(ci  +  C2),  02  = -^ +  Ci +  -C2,  =  4+^1  ■“‘s'- 
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The  error  to  leading  order,  is  given  by  (-3  +  28ci  +  C2)h^/-'*/360.  Note  that  ci  =0,  C2  = 
6/5  yields  the  standard  fourth  order  Fade  scheme  for  the  second  derivative. 

Sixth  order 

Matching  terms  up  to  yields, 

~  ~  33  ~  19 

Co  =  — 6  +  54ci,  C2  =  3  —  28ci,  «2  —  r  — ^ci,  6i  =  — -  +  -ci.  (7c) 

o  Z  o  Z 

The  resulting  error  to  leading  order  is  (1/20160  +  3ci/560)/i® //'*”  .  Note  that  Ci  =  3/44 
yields  the  standard  sixth  order  Fade  scheme  for  the  second  derivative. 

Eighth  order 

Matching  terms  up  to  /"***  yields, 

~  _  13  ~  _  1  ..  88  ..23  ~  1 

^0  2  ’  108’  27’  18’  ^  6* 

The  resulting  leading  order  error  is —h®/f/453600. 

2.3:  The  scheme 

The  interior  scheme  involves  applying  the  equations  derived  in  sections  2.1  and  2.2  at 
each  node.  The  resulting  system  of  2N  equations  is  then  solved  to  obtain  //  and  /".  Of 
the  various  schemes  obtained,  two  schemes  are  discussed  in  detail  below.  These  eure  the 
sixth  order  scheme  with  ci  =  ci  =  0,  and  the  eighth  order  schemes.  These  schemes  have 
the  same  stencil  width  as  the  standard  fourth  and  sixth  order  Fade  schemes.  A  detailed 
comparison  between  these  schemes  and  the  standard  Fade  schemes  is  therefore  performed. 
The  appendix  presents  the  schemes  in  matrix  form,  for  completeness. 

Sixth  order  C-D  scheme  (c^  =  Ci  =01 

rfi-,  + 16/; + 7/;+, + ft(/;i,  -/;;,)  =  ^(/i+i  -  /i-o-  (sa) 

®(/;+l  “  fi-l)  ~  ~  ^fi  +  fi+l)  =  —  2fi  +  /i+i).  (Sb) 

Eighth  order  C-D  scheme 

8i/;_.  + 108/; +51/;+,+ 9/.(/;i.  -/;;,)  =  ^(/i+,  -  /,-. )  -  (9„) 

i3«(/;+.-/;-,)-Mi8/;'..-io8/;'+i8/;;,)  =  -^s±3+Ai+?p(/,+,+/,_,)_^;.. 

(95) 

Standard  fourth  order  Fade 
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fr-i  +  iofi+JUi  =  ^Ui-i-zfi+fi+i)-  (lot) 

Standard  sixth  order  Fade 

fU  +  Vi  +  =  ^(/-»  -  /■-! )  + 

2/!l.+ll/r  +  2/;;,  =  p(/i-,-2/i+/i+,)+j|j(/i-2-2/i+/i+2)-  (lit) 

The  expressions  for  the  first  and  second  derivative  are  seen  to  be  independent  in  the 
standard  Fade  schemes  (equations  10  and  11).  Obtaining  the  first  and  second  deriva¬ 
tives  using  the  standard  Fade  schemes  therefore  involves  separately  inverting  two  tridi¬ 
agonal  matrices  with  band  length  of  N.  By  comparison,  the  first  and  second  deriva¬ 
tives  are  coupled  in  the  C-D  schemes.  The  vector  of  unknowns  is  therefore  of  length 
2JV;  Note  that  for  the  same  stencil  width  as  the  Fade 

schemes,  the  C-D  schemes  are  two  orders  higher  in  accuracy.  This  is  achieved  at  the 
cost  of  inverting  a  matrix  that  has  seven  bands  instead  of  three.  However,  although  the 
band  width  is  increased  from  three  to  seven,  the  inversion  yields  both  the  first  and  second 
derivatives.  A  more  systematic  cost  comparison  with  the  Fade  schemes  is  performed  in 
section  6. 

3.  FOURIER  ANALYSIS  OF  THE  DIFFERENCING  ERROR 

Fomrier  analysis,  and  the  notion  of  the  ‘modified  wavenumber’  provides  a  convenient 
means  of  quantifying  the  error  associated  with  differencing  schemes  [2].  Consider  the  test 
function  fj  =  on  a  periodic  domain.  Discretize  the  function  on  a  domain  of  length 
27r,  using  a  uniform  mesh  of  N  points.  The  mesh  spacing  is  therefore  given  by  h  =  25r fN. 
The  exact  values  of  the  first  and  second  derivative  of  /  axe  and  — .  However, 

the  munerically  computed  derivatives  will  be  of  the  form,  and  —  A;"^e**®-' .  The 

variables  k'  and  are  fimctions  of  k  and  h,  and  are  called  the  modified  wavenumber 
for  the  first  and  second  derivative  operator  respectively.  The  difference  between  k'  and 
Jfe,  and  and  provides  the  differencing  error.  The  modified  wavenumbers  for  the 
coupled-derivative  schemes  are  derived  and  compared  to  the  standard  Fade  schemes,  in 
sections  3.1  and  3.2. 

S.l:  Modified  wavenumber  for  the  standard  Fade  schemes 

The  modified  wavenumbers  for  the  standard  Fade  schemes  are  given  by  Lele  [1]  as 
follows: 

First  derivative 

asinkh +  b/2sin.2kh  \ 

kh  =  - :: — — n -  (12a) 

1  -1-  2a  cos  kh 

where  a  =  1/4,  a  =  3/2,  and  6  =  0  for  the  fourth  order  Fade  scheme.  For  the  sixth  order 
Fade  scheme,  a  =  l/3,a  =  14/9  and  6  =  1/9. 
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Second  derivative 


^,,2^2  _  2a(l  —  coskh)  +  b/2{l  —  cos2kh) 
1  +  2a  cos  kh 


(126) 


where  a  =  1/10, a  =  6/5  and  6  =  0  for  the  fourth  order  Fade  scheme.  For  the  sixth  order 
Fade  scheme,  a  =  2/11,  a  =  12/11,  £ind  6  =  3/11. 


S.2:  Modified  wavenumber  for  the  C-D  schemes 


The  modified  wavenumbers  for  the  C-D  schemes  are  given  below.  As  seen  in  sections 
2.1  and  2.2,  the  sixth  and  eighth  order  schemes  are  members  of  the  following  two-equation 
family  of  schemes: 

/; + “AsWi +/;-i) + hhuui  - f'U) = - /<-.) + - f,-2h  (iso) 


®2(/:+l  — /i-l)+6(6i^:_l  +  +/i-l)  +  -^(/t+2+/i-2)-  (136) 

The  constants  in  the  above  equations  axe  zis  follows: 

Sixth  order  scheme 

C4  =0,  Cl  =  7/16,  62 

a 

ci  =  0,  C2  =  3,  ^  =  —6, 

Eighth  order  scheme 

C4  =  -1/108,  ax  =  17/36,  62  =  -1/12,  C3  =  107/108 

=  -1/108,  C2  =  88/27,  cq  =  -13/2,  2?2  =  23/18,  61  =  -1/6.  (15) 

Equations  13a  and  136  are  used  to  obtain  the  modified  wavenumbers  as  follows.  Con¬ 
sider  the  function,  fi  =  e**^’  on  a  periodic  domain.  Using  the  relations,  =  //e^*** 
and  =  /"e^***,  equations  13a  and  136  become, 

//(I  -h  2ai  cos  kh)  +  ff{i2hb2  sinkh)  =  i-^{cz  sinkh  -f  C4  sin  2A;h).  (16a) 

f-  {i2d2  sin  kh)  +  f-'{h  +  2hbi  cos  kh)  —  ^{co  +  2c2  cos  kh  •+•  2ci  cos  2kh).  (166) 

Equations  16a  and  166  may  be  solved  for  //  and  f^'.  The  resulting  expressions  are  of  the 
form  ik'fi  and  —k"'^fi  where  the  modified  wavenumbers  (after  some  rearrangement)  are 
given  by  the  following  expressions: 


=  -1/16,  C3  =  15/16 

02  =  9/8,  61  =  -1/8.  (14) 
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FIGURE  1  :  The  modified  wavenumber  for  the  first  derivative.  The  C-D  schemes  are  compared  to 

the  standard  Fade  schemes.  - (Exact), - -  (C-D:  eighth  order), . .  (C-D:  sixth  order), 

- (Sixth  order  Fade), - (Fourth  order  Fade). 


,  C3  -h  2C461  -  C062  +  2(0361  -1-  C4  -  62C2)  cos  kh  +  2(0461  -  6201)  cos  2kh 

k  h  =  2smkh - ~ ; - 7- - ~ . - 777 

1  -f-  2c(i6i  -h  20262  "h  2(61  -}"  fli)  cos  kh  -f-  2(ai6i  0262)  cos  2kh 


k"\‘^  =  - 


Oq  -f-  2ctiC2  -f-  202^3  -|-  2(C2  "h  OiCq  202^4)  COS  kh 

1  -|-  2ci6i  -|-  20262  "h  2(61  -f-  oi)cos  kh  -|-  2(oi6i  (i2b2)cos2kh 
2(ci  -t-  0102  —  a2C3)  cos  2kh  -f  4(0101  —  Q2C4)  cos  kh  cos  2kh 
l-|-2oi6i  -1-20262  -f2(6i  +oi)cosfch-f  2(oi6i  —a2h2)cos2kh 


(17a) 


(176) 


3.3:  Evaluation  of  first  derivative 

The  modified  wavenumbers  for  the  first  derivative  are  shown  in  figure  1.  The  C-D 
schemes  are  seen  to  follow  the  exact  solution  more  closely  than  the  standard  Fade  schemes. 
Recall  that  the  sixth  order  C-D  scheme  has  the  same  stencil  width  as  the  fourth  order  Fade, 
while  the  eighth  order  C-D  scheme  has  the  same  stencil  width  as  the  sixth  order  Feide.  In 
spite  of  its  smaller  stencil,  the  sixth  order  C-D  scheme  is  seen  to  have  lower  error  than  the 
sixth  order  Fade.  A  more  quantitative  comparison  of  the  schemes  is  provided  in  table  3. 
The  fractional  error  in  the  first  derivative  may  be  defined  as 


\k'h  -  fch| 


(18) 
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10° 

10 

lO'^ 

10-« 

10-8 


Points  per  wave 

FIGURE  2:  The  percentage  error  in  the  first  derivative  as  a  function  of  the  resolution.  The  C-D 

schemes  are  compared  to  the  standard  Fade  schemes. - -  (C-D:  eighth  order), - (C-D:  sixth 

order), .  (Sixth  order  Fade),  — - (Fourth  order  Fade). 


Points  per  wave 

FIGURE  3:  The  ratio  of  the  error  in  the  first  derivative  between  the  C-D  schemes  and  the  standard 

Fade  schemes  as  a  function  of  the  resolution.  - (C-D  8  /  Fade  6), - (C-D  6  /  Fade  4), 

.  (C-D  6  /  Fade  6). 
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e  =  0.1 

^^1 

Fade  4 

0.59 

0.35 

0.20 

Fade  6 

0.70 

0.50 

0.35 

C-D  6 

0.75 

0.58 

0.42 

C-D  8 

0.81 

0.66 

0.53 

TABLE  3:  A  comparison  of  the  resolving  efficiency  of  the  C-D  schemes  to  the  Fade  schemes. 


iV  =  4 

N  =  8 

Fade  4 

4.51  % 

2.3  X  10"^  % 

Fade  6 

0.97  % 

1.2  X  10-2  % 

C-D  6 

0.36  % 

3.1  X  10-2  % 

C-D  8 

0.06  % 

1.1  X  10-**  % 

TABLE  4:  The  percentage  error  in  the  first  derivative,  as  a  function  of  the  number  of  points  per 
wave  {N),  The  C-D  schemes  are  compared  to  the  standard  Fade  schemes. 

Figure  1  shows  that  the  error  increases  as  kh  increases.  A  measure  of  the  accuracy  or 
‘resolving  ability’  of  the  schemes  is  therefore  provided,  by  specifying  a  maximiun  value  for  e, 
and  estimating  the  fraction  of  the  entire  range  of  wavenumbers  for  which  this  requirement 
is  met.  This  quantity  is  termed  ‘resolving  efficiency’  by  Lele  [1],  and  is  a  hmction  of  the 
specified  tolerance  on  the  error.  Table  3  compares  the  resolving  efficiency  of  the  C-D 
schemes  to  the  standard  Fade  schemes.  The  C-D  schemes  are  seen  to  be  noticeably  more 
accurate.  In  fact,  of  the  different  compact  schemes  considered  by  Lele,  the  only  scheme 
that  outperforms  the  eighth  order  C-D  scheme  is  the  pentadiagonal  tenth  order  scheme 
(designated  ‘i’  by  Lele).  The  pentadiagonal  scheme,  however,  has  a  stencil  of  five  points 
on  the  left  hand  side,  and  7  on  the  right. 

The  modified  waveniunber  may  be  used  to  determine  the  error  as  a  function  of  the 
resolution.  Consider  the  case  where  fc  =  1;  i.e.,  we  have  one  wave  of  wavelength  A  =  2w. 
The  mesh  spacing,  h  is  given  by  h  =  2'k/N  =  X/N.  kh  is  therefore  equal  to  X/N,  the 
reciprocal  of  the  number  of  points  per  wavelength.  The  percentage  error  in  the  first 
derivative  may  be  computed  as  a  function  of  the  resolution,  using  kh  —  2Tr fN,  and  error= 
lOOlfc'h  —  kh\/kh.  Figure  2  compares  the  C-D  schemes  to  the  standard  Fade  schemes. 
Note  that  all  the  schemes  show  100%  error  for  the  two-delta  waves  (two  points  per  wave). 
This  is  because  the  symmetry  of  the  schemes  forces  k'h  to  zero  for  two-delta  waves.  The 
C-D  schemes  are  seen  to  have  noticeably  smaller  error  than  the  standard  Fade  schemes. 
Further  indication  of  this  is  provided  in  figure  3,  where  the  ratio  of  the  error  between  the 
C-D  schemes  and  the  Fade  schemes  is  shown.  Table  4  docmnehts  the  percentage  error  in 
the  first  derivative,  for  resolutions  of  4  and  8  points  per  wave.  The  C-D  schemes  are  seen 
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kh 


FIGURE  4:  The  modified  wavenumber  for  the  second  derivative.  The  C-D  schemes  are  compared  to 

the  standard  Fade  schemes.  - (Exact), - (C-D:  eighth  order), .  (C-D:  sixth  order), 

- (Sixth  order  Fade), - (Fourth  order  Fade). 

to  represent  even  fotur  delta  waves  with  an  accuracy  of  0.4%  and  0.06%,  respectively. 

S.4:  Evaluation  of  second  derivative 

Modified  wavenumbers  for  the  second  derivative  are  shown  in  figure  4.  The  C-D 
schemes  are  seen  to  be  noticeably  more  accurate  at  the  higher  wavenumbers.  Interestingly, 
for  the  C-D  schemes  is  greater  than  the  exact  solution  for  certmn  wavenumbers. 
This  is  in  contrast  to  the  standard  Fade  schemes,  whose  modified  wavenumber  is  eilways 
less  than  the  exact  solution.  However,  this  aspect  of  the  C-D  schemes  does  not  impact  the 
accuracy.  As  shown  in  figures  5  and  6,  the  C-D  schemes  are  more  accurate  than  the  Fade 
schemes  with  the  same  stencil  width. 


Similar  to  the  first  derivative,  a  resolving  efl&ciency  may  be  defined  for  the  second 
derivative,  as  the  fraction  of  the  wavenumber  range  for  which  the  error. 


e  = 


(19) 


is  less  than  a  specified  tolerance.  The  resolving  efficiency  is  tabulated  in  table  5.  Note  that 
the  requirement  that  e  be  less  than  0.1,  is  met  over  the  entire  range  of  waveniunbers.  The 
second  derivative  computed  using  the  sixth  order  C-D  scheme  is  slightly  more  accurate 
than  the  sixth  order  Fade  scheme,  while  the  eighth  order  C-D  scheme  is  noticeably  more 
accurate  than  the  standard  Fade  schemes.  Table  6  shows  the  percentage  error  in  the 
second  derivative,  as  a  function  of  the  resolution.  As  was  observed  for  the  first  derivative, 
the  sixth  and  eighth  order  C-D  schemes  represent  even  four-delta  waves,  to  an  accuracy  of 
about  0.4%  and  0.1%  respectively. 


# 


# 


m 


m 
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FIGURE  5:  The  percentage  error  in  the  second  derivative  zus  a  function  of  the  resolution.  The  C-D 

schemes  are  compared  to  the  standard  Fade  schemes. - (C-D:  eighth  order),  (C-D:  sixth 

order), .  (Sixth  order  Fade), - —  (Fourth  order  Fade). 


FIGURE  6:  The  ratio  of  the  error  in  the  second  derivative  between  the  C-D  schemes  and  the  standard 

Fade  schemes  as  a  function  of  the  resolution.  - -  (C-D  8  /  Fade  6),  (C-D  6  /  Fade  6), 

.  (C-D  6  /  Fade  4). 
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e  =  0.1 

e  =  0.01 

€  =  0.001 

Fade  4 

0.68 

0.39 

0.22 

Pad4  6 

0.80 

0.55 

0.38 

C-D  6 

1.00 

0.57 

0.39 

C-D  8 

1.00 

0.67 

0.50 

TABLE  5:  Comparison  of  resolving  efficiency  of  the  C-D  schemes  to  the  Fade  schemes. 


II 

II 

00 

Fade  4 

2.73  % 

1.6  X  10-1  % 

Fade  6 

0.52  % 

7.41  X  10"®  % 

C-D  6 

0.44  % 

6.16  X  10-3  % 

C-D  8 

0.09  % 

2.84  X  10-^  % 

TABLE  6:  The  percentage  error  in  the  second  derivative,  as  a  function  of  the  number  of  points  per 
wave  (iV).  The  C-D  schemes  are  compared  to  the  standard  Fade  schemes. 


# 


4.  STABILITY  LIMITS  OF  INTERIOR  SCHEME 

This  section  outlines  the  restrictions  imposed  by  Cauchy  stability  on  the  time  step, 
when  the  C-D  schemes  axe  used  with  Runge-Kutta  time  advancement.  The  model  wave 
and  diflEusion  equations  are  solved, 
periodic  dom2iin: 


The  above  equation  is  solved  by  the  method  of  lines.  Let  u  =  ue**^*.  Spatial  discretization 
leads  to  a  set  of  ODEs  of  the  form: 

^  k'hu.  (21) 

at  h 

The  above  equation  is  of  the  form  dy/dt  =  Xy.  It  is  easily  shown  that  numerical  stability 
requires  that 

h  -  ^  ^ 

where  {k'h)jaa.x  denotes  the  maximum  value  of  the  modified  wavenumber  for  the  first 
derivative,  and  (AiA<)inax  denotes  the  upper  bound  imposed  by  numerical  stability,  when 
the  ODE  dy/dt  =  iXiy  is  numerically  integrated.  (AiA<)inax  has  values  of  0,  \/3  and 
2.85  when  the  standard  second,  third  and  fourth  order  Runge-Kutta  schemes  [19]  are 


Consider  the  one-dimensional  wave  equation  on  a 
du  du 
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RK2 

RKZ 

RK4 

Fade  4 

0 

1.0 

1.645 

Faded 

0 

0.871 

1.433 

C-D  6 

0 

0.815 

1.341 

C-D  8 

0 

0.759 

1.249 

Fourier 

0 

0.551 

0.907 

TABLE  7:  The  maximum  CFL  number  allowed  by  numerical  stability. 


RK2 

RKZ 

RKA 

Fade  4 

0.333 

0.417 

0.483 

Fade  6 

0.292 

0.365 

0.423 

C-D  6 

0.208 

0.260 

0.302 

C-D  8 

0.205 

0.256 

0.297 

Fomier 

0.203 

0.253 

0.294 

TABLE  8:  The  maximum  uAtfh?  allowed  by  numerical  stability. 


used  for  time  advancement.  Table  7  lists  the  corresponding  bounds  on  the  CFL  number. 
As  expected,  the  improved  accuracy  at  the  higher  wavenumbers,  reduces  the  maximum 
allowable  CFL  number. 


Similarly,  upper  bounds  on  can  be  obtained  when  the  one-dimensional  diflEu- 


sion  equation. 


du  _  d^u 
dt  ^  dx^ 


(23) 


is  numerically  solved  on  a  periodic  spatial  domain.  Table  8  lists  the  obtained  bounds, 
when  the  C-D  schemes  are  used  with  Rtmge-Kutta  time  advancement.  The  accuracy  of 
the  C-D  schemes  for  the  two-delta  waves  {kh  =  tt)  results  in  the  viscous  restriction  on  the 
time  step  being  nearly  the  same  as  that  for  a  Fourier  spectral  method. 


6.  BOUNDARY  SCHEMES 

Consider  a  spatial  domain  that  is  discretized  by  using  N  points  (including  those  at 
the  boimdaries).  Equations  8  and  9  show  that  the  sixth  order  scheme  can  be  applied  from 
j  =  2  to  AT  -  1,  while  the  eighth  order  scheme  can  be  applied  from  j  =  3  to  iV  -  2. 
For  problems  with  periodic  boundary  conditions,  the  periodicity  of  the  solution  may  be 
used  to  apply  the  same  equations  at  the  boundary  nodes  (see  the  appendix).  However,  for 
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LHS 

RHS 

fl 

0 

Cl  +  C2  -1-  C3  -h  C4 

fi 

do  ^ 

C2  +  2c3  +  3C4 

fi' 

h[ai  -j-  60  +  &i) 

h{c2+2^  C3-f-32  c4)/2! 

fi" 

h^{ai/2\-\-bi) 

(c2  -1-2®  C3  3®  C4)  /3! 

riv 

h^{ail3\  +  bil2l) 

/i®(c2-1-  2^  C3  +  3'‘  C4)/4! 

fi: 

h^{aj4\  +  bi/Z\) 

/i^(c2-H2®  C3-1-3®  C4)/5! 

fV 

(ai/5!  -4-  fei/4!) 

h®  (c2  -f-  2®  C3  -{-  3®  C4)  /6! 

TABLE  9:  Taylor  table  obtained  for  the  boundary  schemes. 


non-periodic  problems,  additional  expressions  are  needed  at  the  boundary  nodes  to  close 
the  system. 

Consider  y  =  1.  The  following  general  expression  may  be  written  for  and 

oo/i  +  fli/a  +  h{bofi  -h  bif")  =  -{■  C2/2  +  C3/3  -h  C4/4).  (24) 

The  corresponding  equation  at  j  =  N  would  be  given  by: 

Oq/n  +  Oi/n-1  ~  +  C2/N-I  +  CifN—2  +  C4/JV-3).  (25) 

The  width  of  the  stencil  on  the  left  hand  side  of  the  above  equation  is  restricted  to 
two.  This  ensures  that  the  number  of  bands  in  the  left-hand  side  matrix  is  still  seven. 
As  was  done  for  the  interior  scheme,  the  constants  in  equation  24  may  be  obtained  by 
expanding  the  terms  in  a  Taylor’s  series,  and  matching  expressions  of  the  same  order. 
Recall  that  we  need  two  independent  equations  at  each  node.  For  the  interior  schemes, 
we  saw  that  60  was  equal  to  0  if  cq  was  equal  to  1,  and  vice-versa.  This  yielded  the  two 
independent  equations.  As  seen  from  tables  1  and  2,  this  relationship  between  oq  and 
60  for  the  interior  schemes  is  a  natural  consequence  of  their  symmetry.  However  for  the 
boundary  schemes,  it  turns  out  that  setting  ao  to  1  does  not  imply  that  60  is  zero,  and 
vice-versa.  The  equation  obtained  when  oq  =  1,  is  the  same  as  that  obtained  when  60  =  1- 
The  following  procedure  is  therefore  used  to  obtain  two  independent  equations.  When 
matching  the  terms  in  the  Taylor  table,  {ao,bo)  is  first  set  equal  to  (1,0).  This  yields  the 
first  equation.  Next,  (ao,6o)  is  set  equal  to  (0, 1).  This  yields  the  second  equation.  These 
expressions  are  derived  below. 

Taylor’s  series  expansion  of  equation  24  yields  table  9.  There  axe  eight  undetermined 
constants  in  equation  1.  Note  that  if  either  of  the  constraints  (ao,6o)  =  (1,0)  or  (0,1) 
is  imposed,  and  the  terms  in  the  Taylor  table  matched,  the  maximum  order  that  can  be 
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obtained  is  five.  The  following  family  of  schemes  is  obtained  by  matching  the  terms  in 
Table  9  to  different  orders.  Consider  first  the  case  where  a©  =  1,  and  bo  —0.  The  resulting 
expressions  may  be  considered  expressions  for  the  first  derivative. 

5.1:  First  equation  (oq  =  l,bo  =  0) 


Expressions  for  the  coefficients  and  the  corresponding  orders  are  given  below. 

Third  order 

Matching  terms  up  to  /•"  yields, 

Cl  = -3  +  C3  +  8c4,  c2  =  3  -  2c3  -  9c4,  ai=2-6c4,  6i  = -- +C3 +  604.  (26a) 


The  leading  order  error  is  then  (1/24  +  C3/I2  +  .  Note  that  C3  —  1/2,  C4  —  0  yields 

the  standard  one-sided,  third  order  compact  scheme  for  the  first  derivative. 

Fourth  order 

Matdiing  terms  up  to  /?"  yields, 

7  1 

Cl  =  - 4c4,  C2  =4-1-15c4,  C3  =  — —  —  i2c4,  ai  =2— 604,  61  =  — 1— 604.  (266) 

2  ^ 


The  error  to  leading  order  is  given  by  (-1/60  -I-  C4/5)h^/i^.  Note  that  C4  =  -1/6  yields 
the  steindard  one-sided,  fourth  order  compact  scheme  for  the  first  derivative. 

Fifth  order 


Matching  terms  up  to  yields. 


The  error  to  leading  order  is  equzil  to  h^fi^fl20. 


(26c) 


5.2:  Second  equation  (ao  =  0,bo  =  1) 

Similar  expressions  are  obtained  when  oq  =  0,  and  60  =  1.  These  expressions  may 
be  considered  relations  for  the  second  derivative.  The  order  of  the  expressions  will  again 
range  from  three  to  five.  However,  due  to  the  second  derivatives  being  multiplied  by  h, 
the  corresponding  order  of  the  second  derivatives  ranges  from  two  to  four.  The  values  of 
the  constants  are  given  below. 

Second  order 


Matching  terms  up  to  yields. 

Cl  =  6 -1- C3 -1- 8c4,  C2  = -6-2c3  -9c4,  ai=-6-6c4,  61  =  2  4- C3 -1- 6C4.  (27a) 


The  error  to  leading  order  is  given  by  (-1/44-  C3/I2  4-  C4)h^/i’'. 
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Third  order 

Matdiing  terms  up  to  /?*'  yields, 

Cl  =  9— 4c4,  C2  =  — 12  +  15c4,  C3  =  3  — 12c4,  ai  =  —6  —  604,  61  =  5— 6C4.  (276) 

The  error  to  leading  order  is  equal  to  (7/60  -f-  C4/5)6^/i .  Note  that  C4  =  —1  yields  the 
standard  one-sided  third  order  compact  scheme  for  the  second  derivative. 

Fourth  order 

Matching  terms  up  to  yields. 

Cl  =  34/3,  C2  =  -83/4,  C3  =  10,  C4  =  -7/12,  ai  =  -5/2,  61  =  17/2.  (27c) 

The  resulting  leading  order  error  is  given  by  — 236^/i'*/60. 

5.3:  Numerical  stability 

The  interior  schemes  outlined  in  section  2.3  are  combined  with  the  boundary  schemes 
of  sections  4.1  and  4.2,  to  close  the  system  of  equations  for  the  first  and  second  derivatives. 
Note  that  the  sixth  order  interior  scheme  may  be  applied  from  j  =  2  to  j  =  N  —  1  zmd 
therefore  only  needs  the  boundary  expressions  at  j  =  1  and  N.  The  eighth  order  interior 
scheme  uses  a  five  point  stencil  on  the  right-hand  side.  It  therefore  can  only  be  applied 
from  y  =  3  to  iV  —  2.  In  this  paper,  if  the  eighth  order  scheme  is  used  in  the  interior,  the 
system  of  equations  is  closed  by  applying  the  sixth  order  scheme  at  j  =  2  and  N  —  1,  and 
the  boimdary  expressions  at  /  =  1  and  N.  These  expressions  are  derived  below. 

Note  that  the  formal  order  of  accuracy  of  the  bovmdary  schemes  is  less  than  the 
interior.  This  is  due  to  the  negative  influence  of  high  order  (wide  stencil)  boundary  closures 
on  the  stability  of  the  overall  scheme.  Past  work  has  shown  that  high  order  boundary 
closiures  can  result  in  munerical  instability  in  hyperbolic  problems.  For  example.  Carpenter 
et  al.  [20]  compute  solutions  to  the  one-dimensional  wave  equation,  and  show  that  the 
standard  foiudih  order  Fade  scheme  (equation  10a)  is  asymptotically  unstable  when  the 
one-sided  foinrth  order  expression, 

/,'+3/J=  ^(-17/1+9/2 +9/3 -A)  (28) 

is  used  at  the  bo\mdary  nodes.  The  third  order  boimdary  expression, 

/i'+2/^  =  ^(-5/i-f-4/2  +  /3)  (29) 


# 


is  shown  to  be  stable.  # 

The  combination  of  the  boundary  and  interior  schemes  is  numerically  tested  for  hy¬ 
perbolic  stability  in  this  section.  The  standard  fourth  order  Runge-Kutta  method  is  used 
for  time  advancement.  The  one-dimensional  wave  equation  is  numerically  integrated  to 
long  times,  and  the  solution  is  examined  for  boundedness  (asymptotic  stability).  Also,  the 
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Time 

FIGURE  7:  Illustration  of  the  asymptotic  instability  of  the  sixth  order  C-D  scheme  with  (4, 3)  closure 
at  the  boundaries.  The  lines  correspond  to  a  CFL  number  of  1.33  while  the  symbols  correspond  to 
a  CFL  number  of  0.1. - ,  •  {N  =  26),  — —  ,  x  (iV  =  51),  - - ,  +  (iV  =  101). 


computational  grid  is  refined  while  keeping  the  CFL  number  fixed,  and  convergence  of  the 
solution  established  (Lax  stability).  Details  of  this  evaluation  are  provided  below. 

Consider  the  one-dimensional  wave  equation. 

Equation  30  is  numerically  solved  over  the  domain  -1  <  x  <  1,  subject  to  the  following 
initial  and  boxmdary  conditions, 

«(x,0)  =  sin27rx,  u(— 1,<)  =  sin27r(— 1  —  t).  (31) 

Note  that  the  exact  solution  to  the  above  equation  is  given  by 


Wexact(®?  0  —  sin27r(a;  t). 


A  uniform  mesh  is  used  for  spatial  discretization.  The  number  of  grid  points  (including 
the  boundaries)  is  set  equal  to  26,  51  or  101.  The  solution  is  then  integrated  to  a  time 
t  =  100.  Note  that  the  solution  travels  one  wavelength  in  one  time  unit,  and  travels  the 
length  of  the  domain  in  two  time  units.  The  L2  error  defined  as. 
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FIGURE  8:  Illustration  of  the  asymptotic  stability  of  the  sixth  order  C-D  scheme  with  third  order 
closure  at  the  boundaries.  The  lines  correspond  to  a  CFL  number  of  1.33  while  the  symbols 

correspond  to  a  CFL  number  of  0.1.  — - ,  •  {N  —  26), - ,  x  {N  =  51), - ,  +  {N  = 

101). 

is  then  examined  for  boundedness. 

Several  combinations  of  the  boimdary  closures,  and  the  interior  scheme  were  examined. 
In  the  following  discussion,  the  notation  [a,  b—c—a,  6]  is  used  to  denote  these  combinations, 
c  denotes  the  order  of  the  interior  scheme,  while  a  2ind  h  denote  the  order  of  the  expressions 
for  the  first  and  second  derivative  at  j  =  1  and  N.  For  example,  the  notation  [3, 3—6—3, 3] 
implies  that  the  sixth  order  scheme  (equation  8)  is  used  in  the  interior,  and  the  third  order 
equations  26a  and  276  are  used  at  the  boundary  nodes.  Note  that  if  c  =  8,  it  is  implied 
that  the  eighth  order  scheme  is  applied  from  j  =  Z  to  N  —  2,  and  the  sixth  order  scheme 
is  applied  at  j  =  2  and  N  —  1. 

The  numerical  evaluations  show  that  the  stability  is  essentially  dictated  by  the  first 
derivative  expression  at  the  boundary.  Schemes  involving  fourth  and  fifth  order  expressions 
for  the  first  derivative,  ».e.,  the  schemes  [4,4  — 6  — 4, 4],  [4,3  — 6— 4,3],  [4,2  — 6— 4, 2],  [5,4  — 

6  —  5,4],  [5,3  —  6  —  5,3],  [5,2  —  6  —  5,2]  were  found  to  be  asymptotically  unstable.  Figtue 

7  illustrates  the  observed  instability  when  the  [4,4  —  6  —  4,4]  scheme  i.e,  fourth  order 
boundary  closure  along  with  a  sixth  order  interior  scheme  is  used.  Note  that  the  L2  error 
is  bounded  at  a  CFL  number  of  1.33  (the  upper  limit  for  stability  of  the  interior  scheme; 
see  Table  7).  However,  the  error  is  seen  to  grow  exponentially  at  a  smaller  CFL  number  of 
0.1.  This  behavior  is  similar  to  that  observed  by  Carpenter  et  al.  [20]  when  the  standard 
fourth  order  Fade  scheme  (equation  10a)  is  used  along  with  a  fourth  order  boundary  closure 
(equation  28).  It  is  a  result  of  the  spatial  discretization  yielding  a  positive  eigenvalue  that 
lies  within  the  stability  envelope  of  the  Runge-Kutta  scheme  at  a  CFL  number  of  1.33. 

Similar  tests  showed  that  combinations  of  the  third  order  expression  for  the  first 
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derivative  (equation  26a)  with  second,  third  and  fourth  order  expressions  for  the  second 
derivative  it  i.e.,  the  schemes  [3, 2  -  c  -  3, 2],  [3, 3  -  c  -  3, 3],  [3, 4  -  c  -  3, 4]  were  stable, 
for  the  sixth  and  eighth  order  interior  schemes.  Figure  8  illustrates  the  stabihty  of  the 
[3, 3  —  6  —  3, 3]  scheme. 


5.4:  Eigenvalue  analysis 

Section  5.3  used  nmnerical  solutions  of  the  wave  equation  to  identify  the  boundary 
closures  that  yielded  stable  solutions  at  long  times.  An  eigenvalue  analysis  is  conducted  in 
this  section  to  confirm  that  these  boundary  closures  do  indeed  yield  asymptotically  stable 
solutions.  Consider  the  wave  equation, 

subject  to  the  inflow  boundary  condition,  m(0,  <)  =  ot  .  Discretize  «  on  a  uniform  grid  of  N 
points  (including  the  boundaries).  The  inflow  condition  implies  that  ui(t)  =  0.  Equation 
34  is  therefore  solved  for  Uj,  where  j  varies  from  2  to  N.  Spatial  discretization  yields  a  set 
of  ODEs  of  the  form: 

^  =  (35) 

dt  h 

where  j  and  k  vary  from  2  to  AT.  M  is  a  iV  -  1  by  AT  -  1  matrix,  and  is  defined  such  that 
u'  =  -MjkUk.  The  eigenvalues  of  M  determine  the  asymptotic  stability  of  the  system 
of  ODEs.  The  requirement  that  the  eigenvalues  of  M  have  negative  real  parts  ensures 
asymptotic  stability.  The  matrix  M  is  obtained  as  follows.  First,  the  condition  ui  =0, 
the  boundary  expressions,  and  the  interior  scheme  are  used  to  eliminate  and  «"  from 
the  system  of  equations  for  the  nodal  derivatives.  The  resulting  system  of  equations  is  then 
rearranged  as  follows.  Recall  that  we  use  two  independent  equations  relating  and  Uj 
at  each  node.  It  is  easily  seen  that  these  two  equations  may  be  expressed  in  the  following 
form:  ^ 

A  u'  +  AB  u"  =  iRi  u,  (36a) 

n 

C  u'  +  AD  u"  =  iR2  U.  (366) 

Note  that  the  above  system  of  equations  is  applied  at  the  nodes  j  =  2  to  AT.  u"  may 
be  eliminated  from  the  above  system  of  equations  to  obtain  an  expression  relating  u  to 
u.  Premultiplying  equations  36a  and  366  by  and  D“^  respectively,  and  subtracting 
yields: 

fB-'A-D-'cV'^KB-'Ri-D-’RiV,  (37) 


f  This  simple  inflow  condition  is  adequate  to  determine  the  inherent  stability  of  the  system.  A 
more  general  inflow  condition,  u(0,<)  =  g{t)  would  simply  yield  a  forcing  term  on  the  right-hand 
side  of  equation  35.  The  stability  of  the  system  would  still  be  governed  by  the  eigenvalues  of  M. 
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implying  that 


{' 


B-^A-D 


-1 


(38) 


Comparison  to  the  relation,  u'  =  —MjkUk  yields  the  following  expression  for  M: 

M  =  ^B-^A  -  D-^C^  (b-^Rx  -  D-^Rz^.  (39) 

The  stability  of  the  (3, 2),  (3, 3)  and  (3,4)  boundary  closures  (section  5.2)  was  tested  for 
both  sixth  and  eighth  order  interior  schemes.  The  nmnber  of  points  N  was  set  equal  to  26, 
51  or  101.  The  matrices  A,B,  C,D,  Ri  and  Rz  were  specified,  and  equation  39  was  used 
to  (numerically)  obtain  M.  An  eigenvalue  solver  from  the  IMSL  library  was  then  used  to 
obtain  the  eigenvalues  of  M.  All  three  boxmdary  closures  were  found  to  yield  eigenvalues 
with  negative  real  parts.  Figures  9  Eind  10  illustrate  the  eigenvalues  obtained  when  the 
(3, 3)  closure  was  used  with  the  sixth  and  eighth  order  interior  schemes  respectively. 


5.5:  The  stable  boundary  closures 

The  stable  boundary  closures  are  summarized  below.  The  following  expressions  are 
used  at  y  =  1.  Equation  25  may  be  used  to  obtain  the  corresponding  expressions  at  j  =  N. 

(3.41  boundary  closure 

The  third  order  expression  for  the  first  derivative  is  combined  with  a  fourth  order  expression 


for  the  second  derivative. 

=  (40..) 

-|/2  +  Kf"  +  ^Si)  =  ^(f/l  -  f  A  +  10/3  -  j|/4)  (406) 

(3. 3)  boundary  closure 

The  third  order  expression  for  the  first  derivative  is  combined  with  a  third  order  expression 
for  the  second  derivative. 

/!+2/2-|A'  =  f(A-/i)  (41«) 

-6/5+M/i'  +  8A')  =  |(3/i-4/2  +  /3)  (416) 

(3.21  boundary  closure 

The  third  order  expression  for  the  first  derivative  is  combined  with  a  second  order  expres¬ 
sion  for  the  second  derivative. 

/i'4-2A-|a'  =  |(A-/x)  (42,.) 

-6A  +  M/i' +  2^')  =  f  (A  -  A)  (426) 


The  matrix  form  of  the  schemes  obtained  with  the  (3, 3)  closure  is  provided  in  the 
appendix  for  completeness.  The  corresponding  matrices  for  the  (3,2)  and  (3,4)  closures 
are  easily  obtained  from  equations  40  and  42. 


•  Real  part 

FIGURE  9:  Eigenvalues  obtained  when  the  (3, 3)  closure  is  used  along  with  the  sixth  order  interior 
scheme.  •  (IV  =  26),  X  (iNT  =  51),  +  (iV  =  101). 


Real  part 


FIGURE  10:  Eigenvalues  obtained  when  the  (3, 3)  closure  is  used  along  with  the  eighth  order  interior 
scheme.  •  {N  =  26),  x  (JV  =  51),  +  {N  =  101). 
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RHS 

LU  solve 

Total 

Fade  4:  first  der. 

l-l-l-|-0  =  2 

2  +  2-l-l  =  5 

3  -f-  3  *1"  1  =  7 

Fade  4:  second  der. 

2-h2-l-0  =  4 

2-f-2-l-l  =  5 

44-44-1  =  9 

Fade  6:  first  der. 

2-H3-1-0  =  5 

2 -f- 2 -1-1  =  5 

44-5  4-1  =  10 

Fade  6:  second  der. 

4-l-5  +  0  =  9 

2-»-2-f-l  =  5 

6  4-74-1  =  14 

Fade  4:  both  ders. 

3-|-3-|-0  =  6 

4  -f-  4  -{-  2  =  10 

74-74-2=  16 

Fade  6:  both  ders. 

6  -H  8  -1-  0  =  14 

4  4-  4  -f-  2  =  10 

10  4- 12  4-  2  =  24 

C-D  6 

3  +  3-|-0  =  6 

12  4-4-1-2  =  18 

15-^74-2  =  24 

C-D  8 

3  -f  7  -H  0  =  10 

12  4-44-2  =  18 

154-114-2  =  28 

TABLE  10:  The  operation  count  per  node  to  compute  the  first  and  second  derivative.  The  entries 
are  of  the  form,  ‘number  of  multiplies  +  adds/subtracts  +  divides  =  total’.  The  overall  cost  is 
obtained  by  multiplying  the  entries  by  the  total  number  of  points,  N. 


6.  COST  COMPARISON 

The  computational  cost  of  the  C-D  schemes  is  compared  to  that  of  the  standard  Fade 
schemes,  in  this  section.  The  standard  Fade  schemes  and  the  C-D  schemes  are  both  of  the 
form, 

Af  =  Bf  (43) 

where  f  =  [. and  A  and  B  are  constant  matrices  that  depend  on 
the  scheme.  For  the  standard  Fade  schemes,  the  vector  f  is  of  length  N,  and  is  either 
equal  to  [. . .  /•_j ,  //,  . .  .]^,  or  [. . .  //ij ,  /f ,  j . .  .]^.  Also,  the  matrix  A  is  tridiag¬ 

onal  with  a  band-length  of  N.  For  the  C-D  schemes,  f  is  of  length  2iV,  and  is  equal 
to  [, . .  /[_j ,  /j'Lj ,  /• ,  The  matrix  A  now  has  seven  bands  (see  the  ap¬ 

pendix),  each  of  length  equal  to  2N. 

At  first  glance,  it  might  appear  as  if  the  C-D  schemes  would  be  significantly  more 
expensive.  However,  this  is  not  the  case.  Although  the  matrix  bandwidth,  and  the  solution 
vector  length  of  the  C-D  schemes  is  twice  that  of  the  standard  schemes,  a  single  evaluation 
yields  both  first  and  second  derivatives.  When  the  cost  of  computing  both  derivatives  is 
estimated,  the  C-D  schemes  are  seen  to  inctur  essentially  the  same  cost  as  the  standard 
Fade  schemes.  This  is  illustrated  below. 

In  using  schemes  of  the  form  given  by  equation  43,  the  common  practice  is  to  perform 
LU  decomposition  of  the  matrix  A  only  once,  and  store  the  L  and  U  matrices.  Compu¬ 
tation  of  the  derivatives  therefore  involves  computing  the  right-hand  side  (B  f),  followed 
by  forward  and  back  substitution.  The  operation  count  associated  with  computing  the 
right-hand  side,  and  solving  the  resulting  system  of  equations  is  tabulated  in  Table  10. 
When  the  cost  of  computing  both  derivatives  is  estimated,  the  C-D  schemes  are  seen  to 
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LU  decomposition 

RHS 

LU  solve 

RHS  -+•  LU  solve 

Fade  4 

2490 

305 

468 

773 

Faded 

2480 

328 

470 

798 

C-D  6 

3590 

251 

412 

653 

C-D  8 

3830 

264 

424 

688 

TABLE  11:  The  time  in  microseconds  to  compute  both  first  and  second  derivatives,  on  a  mesh  of 
128  points.  All  computations  use  LAPACK  routines  for  banded  matrices  on  a  CRAY  C90  in  scalar 
mode. 

involve  the  same  number  of  divides, t  and  add/subtracts  as  the  standard  Fade  schemes 
with  the  garnp  stencil  width.  The  only  increase  in  the  number  of  operations  involves  the 
number  of  multiplies:  the  eighth  order  scheme  has  1.5  times  the  ntunber  of  multiplies  as 
the  sixth  order  Fade,  while  the  sixth  order  scheme  has  twice  the  number  of  multiplies  as 
the  fointh  order  Fade.  As  shown  below,  this  increase  in  the  number  of  multiplies  is  not 
very  significant.  A  cost  evaluation  was  performed  on  a  CRAY  C90  in  scalar  mode,  using 
LAFACK  routines  for  the  LU  decomposition,  and  the  solution  of  the  LU  decomposed 
system.  The  LAFACK  routines  took  advantage  of  the  banded  structure  of  the  coefficient 
matrix.  The  function  /  =  sin(x)  was  discretized  using  a  uniform  mesh  of  128  points  on  a 
domziin  of  length  equal  to  2tt.  Individual  routines  computed  the  right-hand  side,  generated 
and  LU  decomposed  the  matrix  A,  and  solved  the  system  of  equations.  Each  of  these  pro- 
cedmres  was  performed  1000  times,  and  the  result  was  averaged  to  determine  the  cost  per 
evaluation.  The  cost  in  microseconds  is  listed  in  Table  11.  The  C-D  schemes  are  seen  to 
incur  essentially  the  same  cost  as  the  standard  Fade  schemes  when  the  cost  of  computing 
both  derivatives  is  considered. 

7.  CONCLUSION 

A  family  of  finite  difference  schemes  for  the  first  and  second  derivatives  of  smooth 
functions  were  derived.  The  schemes  are  Hermitian  and  symmetric,  and  may  be  consid¬ 
ered  an  extension  of  the  standard  Fade  schemes  described  in  [1].  They  axe  different  from 
the  standard  Fade  schemes,  in  that  the  first  and  second  derivatives  are  simultaneously 
evaluated.  Foiirier  zmalysis  was  used  to  compare  the  proposed  schemes  to  the  standard 
Fade  schemes.  For  the  same  stencil  width,  the  proposed  schemes  were  shown  to  be  two 
orders  higher  in  accuracy,  and  have  significantly  better  spectral  representation.  Numeri¬ 
cal  solutions  to  the  one-dimensional  wave  equation,  and  eigenvalue  analysis  were  used  to 
demonstrate  the  numerical  stability  of  the  schemes.  The  computational  cost  of  the  pro¬ 
posed  schemes  was  assessed,  and  the  cost  of  computing  both  derivatives  was  shown  to  be 
essentially  the  same  as  the  standard  Fade  schemes. 

f  The  divides  in  the  LU  solve  may  be  replaced  by  multiplies  if  desired,  by  storing  the  inverse  of 
the  coefficients  of  the  L  and  U  matrices. 
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Considering  that  the  Navier  Stokes  equations  require  both  first  and  second  derivatives 
of  most  flow  variables,  the  proposed  schemes  appear  to  be  attractive  alternatives  to  the 
standard  Fade  schemes  for  computations  of  the  Navier  Stokes  equations. 
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APPENDIX 


The  schemes  are  presented  in  matrix  form  below.  Both  periodic  and  non-periodic 
boundaries  are  considered. 


The  sixth  order  scheme  on  a  periodic  domain  is  given  by: 


Eighth  order  scheme:  periodic 
The  eighth  order  scheme  on  a  periodic  domain  is  given  by: 
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1 

h 


107(/2-/jv)-(/3-/jV-i) 

-(/3+/jv-i)+352(/2+/jv)-702/i 


107(/,+i-/._i)-(/.+2-/.-2) 

-(/.•+2+/.-2)+352(/.+i+/i_i)-702/i 


107(/i-/jv_i)-(/2-/Ar-2) 

— (/2+/n-2)+352(/i+/iv-i)— 702/n 


The  domain  is  non-periodic.  The  eighth  order  interior  scheme  is  used  at  the  nodes 
j  =  Z  to  N  —  2.  The  sixth  order  interior  scheme  is  used  at  j  =  2  and  N  —  I,  and  third 
order  boundary  expressions  (equation  41)  are  used  at  j  =  1  and  N.  The  resulting  scheme 
is  given  by: 
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Hh-fi) 

9/1-12/2+3/3 

15(/3-/i) 

24(/i-2/2+/3) 


107(/.+i-/._i)-(/.+a-/.-2) 

-(/.•+2+/.-2)+352(/.+i+/._x)-702/i 


15(/jv— /w-a) 

24(/Ar-2— 2/iv_i+/iv) 

3(/jv— /jv-i) 

— 9/Ar+12/)v_i— 3/iv-2 

The  expressions  provided  in  section  5.5  may  be  used  to  obtain  the  matrices  corre¬ 
sponding  to  the  (3, 2)  emd  (3, 4)  boundary  closures. 
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Direct  numerical  simulation  and  inviscid  linear  analysis  are  used  to  study  the  in* 
teraction  of  a  normal  shock  waye  with  an  isotropic  turbulent  field  of  vorticily  and 
entropy  fluctuations.  The  role  of  the  upstream  entropy  fluctuations  is  emphasized. 
The  upstream  correlation  between  the  vorticity  and  entropy  fluctuations  is  shown 
to  strongly  influence  the  evolution  of  the  tui^ulence  across  the  shock.  Negative 
upstream  correlation  between  if  and  T'  is  seen  to  enhance  the  amplification  of  the 
turbulence  kinetic  energy,  vorticity  and  thermodynamic  fluctuations  across  the  shock 
wave.  Positive  upstream  correlation  has  a  suppressing  effect  An  explanation  based 
on  the  relative  effects  of  bulk  compression  and  baroclinic  torque  is  proposed,  and  a 
seeing  law  is  derived  for  the  evolution  of  vorticity  fluctuations  across  the  shock.  The 
validity  of  Morkovin’s  hypothesis  across  a  shock  wave  is  examined.  Linear  analysis 
is  used  to  suggest  that  shock-front  oscillation  would  invalidate  the  relation  between 
Urms  and  T^ms.  as  expressed  by  the  hypothesis. 


1.  introduction 

The  .  interaction  of  shock  waves  with  turbulrait  boundary  layers  has  received  con¬ 
siderable  attention  over  the  past  five  decades.  There  have  been  several  experimental 
studies  (Green  1970;  Femholz  &  Hnley  1981;  Settles  &  Dodson  1994;  Dolling  1993) 
of  the  flow  in  a  compression  comer,  normal  shock/boundary  layer  interaiction  aiid 
more  recently,  the  interaction  of  shock  waves  with  three-dimensional  boundary  layers. 
Th^  experiments  testify  to  the  complex  flow  field  associated  with  tl^'  interaction, 
Incidence  of  the  shock  wave  is  shown  to  strongly  affect  both  the  mean  flow,  and 
the  turbulent  fluctuations.  The  Reynolds  stresses  and  temperature  fluctuations  in  the 
boundary  layer  are  seen  to  amplify  across  the  shock  wave.  Significant  unsteadiness 
of  the  shock  wave  is  observed  when  the  boundary  layer  separates.  The  motion  of 
the  shock  wave  seems  to  be  correlated  with  the  upstream  pressure  fluctuations.  Also, 
high  levels  of  wall-  pressure  fluctuations  are  observed  in  the  vicinity  of  the  shock. 

The  importance  of  bulk  compression  in  the  evolution  of  the  turbulence  across  the 
shock  was  no^  by  Bradshaw  (1974),  Dussauge  &  Gaviglio  (1981)  and  DeWeve, 
Gouin  &.  Gaviglio  (1982).  Amsotropy  of  the  upstream  turbulence  was  identified  as 
a  significant  factor  by  Jacquin,  Blin  &  Geffroy  (1991)  and  Mahesh,  Lele  &  Moin 
(1993).  The  role  of  upstream  acoustic  waves  was  studied  by  Mahesh  et  d.  (1995).  By 
comparison,  the  role  of  entropy  fluctuations  in  shock  wave/boundary  layer  interaction 
seems  to  be  under-appreciated.  For  mcample,  Smits  &  Muck  (1987)  suggest  that  the 
primary  mechanism  of  turbulence  amplification  in  compression-comer  flow  is  inviscid 
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bulk  compression.  The  adverse  pressure  gradient,  concave  streamline  curvature,  and 
mean  compression  doymstream  of  the  comer  are  found  to  further  enhance  turbulence 
levels.  Also,  unsteadiness  of  the  shock  front  is  identified  as  becoming  important 
across  strong  shocks. 

No  mention  is  made  of  the  likely  infiuence  that  upstream  entropy  fluctuations 
would  have  on  the  evolution  of  the  turbulence.  Morkovin’s  (1961)  hypothesis, 

^=-^  =  (7 -DM--,  (U) 

suggests  that  in  addition  to  being  correlated,  the  intensity  of  the  vortical  field  and  ea- 
tropy  Actuations  in  a  constant-pressure  boundary  layer  are  of  comparable  magnitude. 
Experimental  data  (Bradshaw  1977)  support  the  above  hypothesis.  Also,  measure¬ 
ments  (Femholz  &  Finley  1981)  show  that  Prms/p  in  a  constant-pressure  boundary 
layer  can  have  values  as  high  as  5%  - 10%.  This  suggests  that  entropy  fluctuations 
rmght  play  an  important  role  in  the  shock  wave/boundaiy  layer  interaction. 

Past  anal^cal  work  (Morkovin  1960;  Chang  1957;  Cuadra  1968)  has  examined 
the  mteraction  of  plane  entropy  waves  with  a  shock.  These  studies  emphasize  the 
pressure  field  that  is  produced  through  the  interaction.  Morkovin’s  one-dimensional 
analyst  s^ws  that  the  pressure  fluctuations  generated  behind  a  detached  shock 
Me  quite  intense,  and  may  affect  the  transition  of  the  boundary  layer  on  the  body. 
Cuadra’s  parametric  study  demonstrates  that  entropy  fluctuations  in  shear  flows  can 
generate  noise  levels  comparable  to  those  generated  by  vortical  fluctuations  Rapid 
distortion  theory  was  used  by  Goldstein  (1979)  to  examine  the  influence  of  upstream 
entropy  fluctuations  on  turbulence  passing  through  a  wind-tunnel  contractioa  His 
Malysis  showed  that  the  entropy  fluctuations  produced  turbulence,  whose  magnitude 
mc^ased  more  rapidly  with  contraction  ratio  than  that  of  the  upstream  imposed 
turbulence.  With  the  exception  of  Hesselink  &  Sturtevant’s  (1988)  experiments,  it 
appears  that  the  interaction  of  a  turbulent  field  of  entropy  fluctuations  with  a  shock 
wave  has  not  been  examined.  Also,  Hesselink  &  Sturtevant  were  interested  in  the 
some  boom  problem.  As  a  result,  their  study  emphasized  the  evolution  of  the  shock 
front,  and  restricted  the  shock  strength  to  Mach  1.1. 

Tks  paper  studies  the  interaction  of  a  normal  shodc  with  a  turbulent  field  of 
wrtidty  and  entropy  fluctuaAns.  Direct  numerical  simulation  (DNS)  and  invisdd 
linear  analysis  are  used  for  this  purpose.  The  focus  is  on  the  influence  of  the  upstream 
en^py  fluctuations.  Shock  waves  of  strength  Mach  U9  and  Mach  1.8  are  computed 
using  DNS,  while  the  linear  analysis  considers  a  range  of  Mach  numbers  from  1  to 
3.  The  pa^rer  is  organized  as  follows.  A  description  of  the  problem  is  prflVdded  in  §2 
Some  details  of  the  DNS  and  linear  analysis  arc  also  summarized.  Section  3 
the  results  and  provides  an  explanation  for  the  rde  of  entropy  fluctuations.  Also,  a 
scaling  law  for  the  evolution  of  vortidty  fluctuations  across  the  shock  is  derived.  ITie 
paper  is  then  concluded  with  a  brief  summary  in  §4. 


2.  Description  of  the  problem 

Direct  numerical  simulation  and  linear  analysis  are  used  to  study  the  interaction 
of  a  normal  shock  wave  with  turbulence.  The  turbulent  velodty  field  upstream  of  the 
shock  \rave  is  isotiopic,  and  the  mean  flow  is  uniform.  Also,  the  upstream  turbulence 
essentally  comprises  a  vortical  velodty  field  and  entropic  thermodynamic  field.  This 
cornpoaAo  ^  exactly  enforced  within  the  linear  analysis,  and  is  only  approximately 
satefied  in  the  DNS.  The  upstream  thermodynamic  field  in  the  DNS  approximately 
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Figure  1.  Schematic  of  the  numerical  simulation. 


satisfies  the  weak  fonn  of  Morkovin’s  hypothesis.  The  analysis  considers  two  different 
spectra  for  the  upstream  thermodynamic  field:  one  that  satisfies  the  strong  form 
of  Morkovin’s  h^othesis  (axisymmetric),  and  the  other  that  is  isotropic.  Note  that 
i/T  is  negative  if  Morkovin’s  hypothesis  holds  upstream  of  the  shock  wave.  For  the 
isotropic  upstream  thermodynamic  field,  the  analysis  considers  two  possibilities:  one 
where  ufT'  is  negative,  and  the  other  where  it  is  positive. 

Section  2.1  outlines  the  relevant  details  of  the  DNS.  The  inflow  and  outflow 
boundary  conditions  are  emphasized.  A  brief  description  of  the  linear  analysis  is  then 
provided  in  §2.2. 


2.1.  Details  of  the  direct  numerical  simulation 

A  schematic  of  the  computed  flow  is  shown  in  figure  1.  Note  that  the  shock  wave  is 
stationary  in  the  mean.  The  governing  equations  are  the  unsteady  three-dimensional 
compressible  Navier-Stokes  equations  in  the  following  conservative,  non-dimensional 
form: 


dp 

d 

dt 

dxi' 

—(oui)  — 

_  d 

dxj 

dEt 

d 

dt  ~ 

dXi 

5x/ 

^  +  — 
dxf  dxj 


(2.1fl) 

(2.1&) 

(2.1c) 


The  viscous  stress  tensor  and  heat  flux  vector  ate  given  by 


(22a) 


ft  dT 
RePrdXi 

The  variable  Et  denotes  the  total  energy,  defined  as  £«  =  p/{y — 1)+ pu(Uill.  Note  that 
the  mean  sound  speed,  density  and  dynamic  viscosity  at  the  inflow  of  the  domain  are 
used  to  non-dimensionalize  velocity,  density  and  viscosity  respectively.  The  reference 
lengthscale  Lr  is  related  to  the  other  reference  variables  by  Re  =  prCrLrl/jir-  The  fluid 
is  assumed  to  be  an  ideal  gas  with  1.4  as  the  ratio  of  specific  heats.  The  dynamic 
viscosity  is  related  to  the  temperature  by  a  power  law  with  0.76  being  the  exponent, 
and  the  Prandtl  number  is  assigned  a  constant  value  of  0.7. 

The  computational  mesh  is  uniform  in  the  directions  transverse  to  the  shodc  wave. 
A  non-uniform  mesh  is  used  in  the  streamwise  direction,  such  that  points  are  clustered 
in  the  vicimty  of  the  shock.  The  following  analytical  mapping  is  used  for  this  purpose. 
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Using  the  variable  s  to  denote  a  uniform  mesh  from  0  to  lx,  the  non-uniform  mesh 
is  given  by 

+  Hs  —  3c/2)  cosh  bc/2' 

X  _  2b  ^  cosh  b(5  — c/2)  cosh  3bc/2 
Lx  ^  CQsh  b(l  -  3c/2)  cosh  bc/2''' 

^  2b  [cosh  b(l- c/2)  cosh  3bc/2 

Values  of  b,r,d  and  c  used  in  the  Mach  1.29  simulations  are  12,  1.95,  0.04  and 
0.281  respectively.  The  Mach  1.8  computation  uses  values  of  12,  1.95,  0.02  and  0.283 
respectively. 

A  combination  of  the  sixth-order  Fade  scheme  (Lele  1992)  and  the  sixth-order  ENO 
scheme  (Shu  &  Osher  1988, 1989)  is  used  to  compute  spatial  derivatives.  The  shock¬ 
capturing  scheme  (ENO)  uses  the  global  Lax-Friedrichs  flux  splitting.  The  degeneracy 
in  the  base  form  of  the  ENO  scheme  is  removed  (Shu  1990)  by  biasing  the  adaptive 
choice  of  stencil  towards  central  differences.  Also,  the  shodc-capturing  scheme  is 
applied  only  in  the  streamwise  (shock-normal)  direction  in  the  vicinity  of  the  shock 
wave.  The  sixth-order  Fade  scheme  is  therefore  used  to  compute  all  spatial  derivatives 
except  the  streamwise  invisdd  fluxes  around  the  shock:  the  ENO  scheme  is  used  to 
compute  those  terms.  Time  advancement  is  performed  using  the  compact-storage 
third-order  Runge-Kutta  scheme  (Wray  1986).  Since  the  turbulence  is  statistically 
homogeneous  tonsverse  to  the  shock  wave,  periodic  boundary  conditions  are  imposed 
in  those  directions.  Turbulent  fluctuations  are  superposed  onto  the  mean  field  at  the 
inflow  boundary  while  non-reflecting  boundary  conditions  are  specified  at  the  exit 
through  use  of  a  ‘sponge’  zone.  Details  of  the  inflow  and  outflow  boundary  conditions 
are  provided  below.  The  Mach  1.29  simulations  are  used  for  illustration;  the  same 
procedure  was  followed  for  the  Mach  1.8  computation. 

Zl.l.  Inflow  turbtdence 

Since  the  flow  upstream  of  the  shock  wave  is  supersonic,  all  the  flow  variables 
are  specified  at  the  inflow  boundary.  Turbulent  fluctuations  in  velocity,  density 
and  pressure  are  superposed  on  the  uniform  mean  flow.  The  turbulence  upstream 
of  the  shock  wave  is  required  to  have  an  isotropic  velocity  field.  When  entropy 
fluc^tions  are  present,  the  upstream  velocity-temperature  correlation  is  required  to 
be  similar  to  that  found  in  adiabatic  boundary  layers,  Le.  is  required  to  be  nearly 
—1.  The  correlation  between  the  velocity  and  temperature  (entropy)  fluctuations  in 
turbulent  boundary  layers  is  a  consequence  of  the  turbulent  velocity  field*stirring  .the 
mean  temperature  (entropy)  gradients  to  produce  temperature  (entropy)  ^uctuations. 
However,  due  to  the  absence  of  a  mean  temptyature  gradient  upstream  of  the 
shock  in  the  computation,  there  is  no  mechanismi  to  sustain  a  negative  velocity- 
temperature  correlation.  This  may  be  anticipated  from  the  independence  of  the 
velocity  and  entropy  fluctuations  in  the  linear  invisdd  limit  for  uniform  mean  flow, 
^e  following  procedure  is  therefore  used  to  generate  the  desired  upstream  turbulence 
in  the  computations.  . 

The  turbulent  fluctuations  are  obtained  from  a  separate  temporal  simulation  of 
isotropic,  drying  turbulence.  The  temporal  simulation  (which'has  periodic  boundary 
conditions  in  all  three  directions)  is  advanced  in  time  until  the  flow  field  is  wefl- 
developed,  ie.  the  velodty  derivative  skewness  5*  =  attains 

(Tavoularis,  Bcmett  &  Corrsin  1978;  Erlebacher  et  al  1992)  a  value  between  -0.4 
and  -^.6.  Typically,  ^  happens  after  a  time  t  ~  X/u^.  Taylor’s  hypothesis  is 
then  invoked,  and  a  single  realization  of  the  developed  field  is  convected  through 
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the  inflow  plane  of  the  spatial  simulation.  Details  of  this  procedure  are  described 
by  Mahesh,  Lele  &  Moin  (1996).  Note  that  Lee,  Lele  &  Moin  (1992)  have  shown 
the  validity  of  Taylor’s  h3rpothesis  in  isotropic  compressible  mrbulence  for  fluctuation 
Mach  number  M,  =  (q^y^/c,  and  turbulence  intensity  as  high  as  0.5  and  0.15 
respectively.  As  shown  in  table  1,  the  turbulence  at  the  inflow  is  within  these  limits. 

The  simulations  with  the  shock  wave  compare  two  different  cases  at  each  Mach 
number.  The  two  cases  differ  in  the  nature  of  the  turbulence  upstream  of  the  shock 
wave.  While  the  upstream  turbulence  in  one  case  (case  A)  is  composed  of  vortical 
fluctuations,  a  combination  of  vorticity  and  entropy  fluctuations  is  present  in  the 
other  (case  B).  The  upstream  entropy  field  in  case  B  is  required  to  approximately 
satisfy  the  weak  form  of  Morkovin’s  hypothesis,  Le.  (1.1)  is  satisfied  in  the  rjn.s.  sense. 
Also,  the  inflow  spectra  and  r.m.s.  levels  of  velocity  in  both  cases  are  required  to 
essentially  be  the  same.  The  same  temporal  simulation  is  therefore  used  to  generate 
the  inflow  turbulence  in  the  two  cases. 

This  is  done  as  follows.  A  temporal  simulation  of  isotropic  turbulence  is  first 
conducted.  The  initial  velocity  field  is  solenoidal  with  energy  spectrum  E{k)  ~ 
and  the  initial  thermodynamic  fluctuations  are  set  to  zero.  The  simulation 
is  advanced  in  time  until  the  velocity  derivative  skewness  attains  its  developed  value. 
An  instantaneous  flow  field  is  then  taken.  The  thermodynamic  fluctuations  in  this 
field  are  nearly  isentropic.  To  generate  a  realistic  turbulent  field  of  vorticity  and 
entropy  fluctuations,  the  following  procedure  is  carried  out  The  pressure  fluctuations 
in  the  stored  flow  field  are  set  to  zero  while  density  fluctuations  that  satisfy  p'/p  = 
(y  -  l)Mfp„i//Uspa,  are  specifiedt  (the  subscript  spat  denotes  spatial).  This  modified 
field  is  then  advanced  in  time.  Statistics  from  the  simulation  are  compared  to  a 
parallel  simulation,  where  the  same  field  without  the  above  modifications  is  advanced 
for  the  same  length  of  time. 

As  expected  from  Kovasznay’s  (1955)  modal  decomposition,  the  entropy  fluctua¬ 
tions  that  are  introduced  do  not  significantly  influence  the  evolution  of  the  velocity 
fluctuations.  After  a  brief  acoustic  transient  (t  MtX/urms),  the  decay  rate  of 
turbulence  kinetic  energy,  and  the  velocity  derivative  skewness  match  those  in  the 
simulation  without  entropy  fluctuations.  However  as  expected,  the  solution  deviates 
from  Morkovin’s  hypothesis  as  it  evolves  in  time.  To  ensure  that  the  weak  formulation 
of  the  hypothesis  is  approximately  satisfied  by  the  turbulence  upstream  of  the  shock, 
an  mstantaneous  realization  is  taken  immediately  after  the  acoustic  transient,  and 
used  to  specify  inflow  turbulence  in  case  B.  A  realization  at  exactly  the  same  instant 
of  time  is  taken  from  the  temporal  simulation  without  entropy  fluctuations  and  used 
to  specify  inflow  turbulence  incase  A.  x 

Figures  2  and  3  present  results  &om  the  temporal  simulation  used  to  genfrate 
inflow  turbulence  in  the  Mach  1.29  simulations.  Note  that  the  solution  at  the  end 
of  the  tem^ral  simulation  is  used  to  generate  the  turbulence  upstream  of  the  shock 
wave.  The  initial  temporal  velocity  field  is  chosen  to  have  fluctuation  Mach  number 
M  =  0.22  Md  microscale  Reynolds  number  Rx  =  Urms^/y  =  39.5.  A 

uniform  mesh  of  81^  points  is  used  on  a  domain  of  length  2n  in  all  tbee  directions. 
The  solution  is  advanced  for  time  t  =  LSSt^,  where  is  a  turbulence  time  scale 
defined  as  the  ratio  of  X  to  Urms  at  f  =  0.  As  shown  in  figure  2,  the  velocity  derivative 
skewness  has  attained  a  value  of  -^.48  by  this  time.  Also,  Mf  and  Rx  have  dropped  to 
0.16  and  22.8  respectively.  Entropy  fluctuations  are  then  introduced,  and  the  solution 

t  If  these  entropy  fluctuations  were  specified  at  r  0,  would  drop  to  very  small  levels  by 
the  time  the  turbulence  is  well^developed. 
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Figum  2.  Temporal  evolution  of  the  velocity  derivative  skevraess  in  the  temporal  decay  of 
isotropic  turbulence; - ,  without  entropy  fluctuations; - ,  with  entropy  fluctuations. 


figure  3.  Tempoi^  evolution  of  r.m.s  values  of  terms  in  Morkovin’s  hypothesis  in  the  temporal 
decay  of  Botropic  turbulence: - .  p^/p;  — -  .  t„/T; . .  [y  -  l]M^wC 


K  advanced  in  time.  At  the  end  of  the  simulation,  the  skewness  has  leattained  its 
deveh>^  value  (figure  2),  and  the  weak  form  of  the  hypothesis  is  atwraximately 
Prm*/?,  '^rm/T  and  (y  -  l)Mj^Urms/Uspat  have  values  of  0.043,  0.043  and 
0.044  restively.  By  comparison,  Pr„s/yp  is  0.01.  Also,  M„  Rx  and'0™„/(»™,  have 
values  of  0.14,  20.6  and  0.085  respectively  (0^,  and  denote  the  rjits.  dilatation 
Md  vomaty  respectively).  The  solution  is  therefore  considered  to  be  dominated 
by  vorticity  and  entropy  fluctuations  that  approximately  satisfy  the  weak  form  of 
Morkovm’s  hypothesis. 


The  solution  at  the  end  of  the  calculations  shown  in  figures  2  and  3  is  used 
to  specify  inflow  turbulence  in  the  Mach  1.29  spatial  simulations.  The  relevant 
parameters  of  the  inflow  turbulence  are  listed  in  table  1.  The  fidelity  of  the  temporal 
smula^ns  was  checked  by  examining  the  energy  spectra  and  two-point  correlations 
of  the  fluctuations.  One-dimensional  spectra  of  the  velocity  field  show  very  good 
agreement  between  the  two  cases,  A  and  B.  The  spectra  show  about  five  decades  of 
drop-ofi,  mdicating  adequate  resolution.  Also,  two-point  correlations  of  the  velocity 
and  density  field  drop  off  to  zero,  indicating  adequate  size  of  computational  domain. 
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2.1.2.  Outflow  boundary  conditions 

Approximately  non-reflecting  boundary  conditions  are  specified  at  the  subsonic 
outflow  boundary.  The  boundary  conditions  comprise  a  ‘sponge  layer’  in  the  stream- 
wise  direction,  followed  by  a  characteristics-based  boundary  condition  (Poinsot  &  Lele 
1992)  at  the  outflow  plane.  Boundary  conditions  involving  a  sponge  layer  have  been 
used  in  the  past  in  a  variety  of  problems:  e.g.  Israeh  &  Orszag  (1981),  Givoli  (1991), 
Colonius,  Moin  &  Lele  (1995).  The  boundary  conditiorrs  with  a  sponge  were  shown  to 
yield  significantly  better  results  than  boundary  conditions  without  the  sponge  layer. 

The  governing  equations  in  the  ‘sponge  layer’  are  modified,  such  that  the  solution  in 
the  layer  is  gently  damped  to  a  reference  solution.  A  term  of  the  form  —o{U—Uref)  is 
added  to  the  right-hand  side  of  governing  equations  over  the  sponge  layer  alone.  £/’„/ 
denotes  the  vector  of  reference  variables  towards  which  the  solution  in  the  sponge 
layer  is  forced.  It  is  obtained  from  the  Ranldne-Hugoniot  equations  for  a  laminar 
shoclc  The  coefficient  a{x)  is  a  polynomial  function;  Le. 

where  Xs  and  L*  denote  the  start  of  the  sponge  and  the  length  of  the  domain 
respectively.  Values  of  A,,  n  and  (L*  —  Xs)/Lx  used  in  our  simulations  are  5,  3  and 
0.14  respectively. 

The  outflow  boundary  conditions  are  evaluated  below.  The  interaction  of  the  Mach 
129  shock  with  a  plane  vorticity-entropy  wave  (45®  incidence  angle)  is  computed.  Lin¬ 
ear  analysis  predicts  that  the  incident  wave  would  generate  downstream-propagating 
acoustic,  vorticily  and  entropy  waves  behind  the  shock.  The  ability  of  the  boundary 
conditions  to  allow  these  waves  to  smoothly  exit  the  domain  is  thus  tested.  The 
mean  flow  parameters,  domain  length  in  the  x-direction,  and  grid  in  the  x-direction 
are  matched  to  those  in  the  turbulent  simulatiotL  The  wavenumber  of  the  upstream 
disturbance  is  set  equal  to  the  wavenumber  at  which  the  energy  spectrum  of  the 
upstream  turbulence  peaks.  The  rms.  level  of  the  upstream  disturbance  is  matched 
to  that  of  the  upstream  turbulence.  The  extent  of  the  domain  in  y  is  set  equal  to  one 
wavelength  of  the  upstream  disturbance,  and  a  grid  of  231  by  16  points  is  used  to 
discretize  the  flow.  The  disturbance  at  the  inflow  boundary  is  given  by  the  real  part  of 
(A  Ic-e)  (see  the  Appendix)  where  the  variables  k,y)i,A„  and  Ae  are  set  equal  to  5, 45% 
0.05  and  0.05  respectively.  The  computation  is  initialized  by  a  numerically  cmuputed  ... 
steady  laminar  shock.  The  disturbance  is  then  introduced  at  the  inflow  boundary. 
Statistics  are  gafliered  over  each  period  of  the  inflow  disturbance,  after  one  domain 
flow-through  time  has  passed.  The  statistics  from  succ^ve  periods  are  compared  to 
diedc  if  initial  transients  persist  After  a  time  Uit/Lx  =  3.5,  the  transient  effects  are 
found  negligible  and  the  statistics  have  converged. 

To  evaluate  the  outflow  boundary  conditions,  the  same  flow  is  computed  on  a 
domain  twice  as  long  behind  the  shock  wave.  1100  points  are  used  in  the  streamwise 
directiorL  The  resulting  resolution  is  greater  than  that  in  the  shorter  domain.  Statistics 
from  the  two  simulations  are  then  compared  to  each  other  and' to  linear  analysis. 
Hgure  4  shows  the  streamwise  evolution  of  the  averaged  kinetic  energy,  vortidty  and 
density  in  the  two  computations.  Ordy  the  ‘useful  region’  (L*  —  x*)  of  the  shorter 
domain  is  showiL  Good  agreement  between  the  two  computations  is  observed, 
indicating  that  the  non-reflecting  nature  of  the  sponge  region  is  acceptable.  Some 
influence  of  the  sponge  region  (maximum  value  about  3.5%)  on  the  statistics  of 
and  is  observed  immediately  upstream  of  the  sponge. 
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Rgure  4.  Spatial  evolution  of  statistics  from  a  computation  of  the  interaction  of  a  Mach  1.29 
shock  wave  with  a  plane  vortoty-entropy  w^.  Results  from  short  and  long  domains  are  both 

rn^cat^  as  follows:  —  - ,  .  ,  - .  (P^m^)/{ql/Vl). 

The  subscnpt  m  denotes  inflow.  ^ 


■Die  spatial  evolution  of  the  statistics  is  not  compared  to  inviscid  linear  analysis 
owing  to  significant  viscous  decay  in  the  computation.  However,  the  amplification  of 
statistics  immediately  across  the  shock  wave  is  compared  to  anal^sJHie  comparison 
shows  that  the  error  (computed  with  respect  to  analysis)  in  and  is 

0.8%,  1.6%,  0.8%  and  0.5%  respectively.  Linear  analysis  shows  that  only  waves 
whose  incidence  angle  is  less  than  a  critical  incidence  an^e  generate  downstream- 
propagating  disturbances  behind  the  shock.  Evanescent  waves  are  generated  at 
hi^er  angles  of  incidence.  In  computing  shock/turbulence  interaction,  the  fraction 
of  incident  waves  that  generate  these  downstream-propagating  disturbances  poses  the 
primary  challenge  to  outflow  boundary  conditions.  The  performance  of  the  outflow 
boundary  for  the  incidence  angle  of  45*  suggests  that  the  boundary  conditions  are 
quite  adequate  for  the  turbulent  simulations. 

In  ev^uating  the  outflow  boimdary  conditions,  it  was  mentioned  that  the  stream- 
wise  gnd  was  the  same  as  that  used  in  the  turbulent  simulation.  This  streamwise 
gnd  is  chosen  as  follows.  The  shock-capturing  scheme  is  used  to  compute  the  cor¬ 
responding  steady  laminar  shock,  and  the  solution  is  evaluated.  The  pNO  schemes 
are  known  (Meadows,  Caughy  &  Casper  1993;  Woodward  &  CoDela  1984;’  Roberts 
1990;  Lindquist  &  Gfles  1991)  to  produce  spurious  entropy  oscillations  when  applied 
to  dowly  moving  shock  waves.  In  contrast  to'^standard  central  difierences,  these 
oscillations  are  bounded;  they  do  not  cause  instability.  The  spurious  oscillations  in 
density  and  temperature  are  typically  the  largest,  those  in  velocity  and  pressure  are 
much  smaller.  To  compute  shock/turbulence  interaction,  we  require  that  the  ampli¬ 
tude  of  these  spurious  oscillations  be  small  compared  to  the  physical  fluctuations. 
The  amplitude  of  these  spurious  oscillations  can  be  reduced  by  increasing  the  shodc 
speed  relative  to  the  mesh  (Meadows  et  al.  1993),  or  by  refining  the  mesh  (for  a 
viscous  computatioii).  The  shock  is  stationary  in  the  mean  in  our  computations. 
The  shock  speed  with  respect  to  the  gnd  is  therefore  determined  entirely  by  the 
upstream  disturbances.  We  therefore  use  the  shock-capturing  scheme  to  compute  the 
corresponding  laminar  shock,  and  refine  the  mesh  until  the  spurious  oscillations  are 
acceptably  small  For  example,  the  mesh  used  in  our  Mach  1 JI9  computations  yields 
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spurious  density  oscUIations  that  are  0.03%  of  the  mean  density  behind  the  Mach  129 
shock  (the  spurious  oscillations  in  temperature,  velocity  and  pressure  are  smaller). 
Table  1  shows  that  the  upstream  intensity  of  density  fluctuations  in  the  turbulent 
simulations  is  ^ater  than  1.2%.  The  above  level  of  spurious  oscillations  generated 
by  the  shock-capturing  scheme  is  therefore  considered  acceptable.  This  streamwise 
grid  is  then  used  to  perform  the  computations. 

21.3.  Simulations  performed 

The  relevant  parameters  of  the  numerical  simulations  are  given  in  table  1.  Note 
that  the  values  at  the  inflow  are  quoted.  Essentially,  the  simulations  consider  the 
interaction  of  low-Reynolds-number  turbulence  with  shock  waves  of  strengths  equal 
to  Mach  1.29  and  Mach  1.8.  Note  that  one  of  the  cases  (free-stream  Mach  number 
2.9,  wedge  angle  8°)  in  the  compression-comer  experiments  of  Smits  &  Muck  (1987) 
yields  a  normal  Mach  number  of  129  if  one  assumes  that  the  entire  flow  is  turned 
across  a  single  shock.  Also,  only  one  simulation  (with  upstream  entropy  fluctuations) 
is  performed  for  the  Mach  1.8  shock.  In  §21  it  was  noted  that  the  shock-capturing 
scheme  was  only  applied  in  the  streamwise  direction  for  the  Mach  129  shock  wave. 
However,  applying  the  shock-capturing  scheme  in  only  the  streamwise  direction  set 
up  unstable  oscillations  for  the  Mach  1.8  shock.  As  a  result,  the  shock-capturing 
scheme  was  applied  in  all  three  directions  in  the  immediate  vicinity  of  the  shock  for 
the  Mach  1.8  shock  wave. 

22.  linear  analysis 

Details  of  the  linear  analysis  are  provided  in  the  Appendix.  Essentially,  the  lin¬ 
earized  Euler  equations  are  used  to  study  the  interaction  between  the  upstream 
turbulence  and  the  shock  wave,  which  is  modelled  as  a  discontinuity.  The  up¬ 
stream  turbulence  (being  homogeneous)  is  represented  as  a  superposition  of  Fourier 
modes  (plane  waves),  each  of  which  independently  interacts  with  the  shock.  The 
analysis  is  therefore  initiated  by  considering  the  interaction  of  a  single  plane  wave 
with  the  shock.  Workers  such  as  Ribner  (1953,  1954)  and  Chang  (1957)  have 
developed  the  procedure  to  obtain  the  disturbance  field  behind  the  shock  and 
the  distortion  of  the  shock  front  Other  references  to  work  using  linear  analy¬ 
sis  are  provided  by  Mahesh  et  al.  (1995).  This  solution  is  then  integrated  over  ’ 
all  the  incident  waves  to  obtain  a  statistical  description  of  the  turbulence  behind 
the  shock  wave.  The  turbulence  statistics  thus  obtained  are  functions  oMhe  "dis-- 
tance-behind  the  shock  wave,  the  Mach  number  of  the  shock,  and  Ihe  incident 
velocity  and  density  spectra.  In  presenting  results -^ipm  the  linear  analysis,  this 
paper  will  occasionally  refer  to  the  variables  Ar  and  Essentially,  A,  teptt- 
sents  the  amplitude  of  the  upstream  density  field,  and  0,  represents  the  sign  of 
the  upstream  velom^— temperature  correlatiotL  If  0^  ~  0,  itT'  is  negative  up¬ 
stream  of  the  shock  wave;  0,  =  180“  implies  that  I/'T'  is  positive  upstream  of  the 
shock. 


3.  Results 

Turbulence  statistics  in  the  simulation  are  computed  by  averaging  over  time  and 
the  y-  and  z-directions.  The  ratios  of  mean  velocity,  temperature  and  pressure 
across  the  shock  wave  are  very  nearly  equal  to  their  corresponding  laminar  values 
(Mahesh  et  al.  1996).  A  small  overshoot  in  pressure  and  temperature  is  observed 
inunediately  downstream  of  the  shock.  The  mean  velocity  exhibits  a  corresponding 
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Case  1.29A 

Case  1.29B 

Case  1.8 

Ml 

1.29 

1.29 

1.8 

Ri 

19.1 

19.1 

19.5 

M, 

0.14 

0.14 

0.18 

(9^)‘^Vl/i 

0.11 

0.11 

0.10 

PrmsrPl 

0.012 

0.042 

0.076 

Wr, 

0.0042 

0.041 

0.076 

Prm/yPx 

0.011 

0.010 

0.018 

n’  1‘^rms  'Rrms 

-0.06 

-0.84 

-0.85 

(231,  81,  81) 

(231,81,81)  1 

(231,  81,  81) 

Table  1.  Parameters  in  the  simulations  performed. 


undershoot,  following  which  it  remains  practically  constant  However,  mean  pressure 
and  temperature  exhibit  a  small  positive  gradient  behind  the  shock  wave.  The 
presence  of  upstream  entropy  fluctuations  has  no  noticeable  effect  on  the  mean  flow 
behind  the  shock.  Entropy  fluctuations  do  however  influence  the  of  the 

shock  wave  as  inferred  from  the  mean  flow  profiles.  The  ‘mean  shock  thickness’  is 
larger  in  the  presence  of  entropy  fluctuations^  Le.  the  mean  gradients  in  the  vicinity 
of  the  shock  ^  smaller.  Given  that  the  ‘mean  shock  thickness’  reflects  the  amplitude 
of  shock  oscillation,  this  indicates  an  increase  in  shock  motion  when  the  upstream 
fluctuations  satisfy  Morkovin’s  hypothesis. 

The  following  sections  present  results  for  the  evolution  of  the  turbulence  across 
the  shock.  The  results  from  the  DNS  and  linear  analysis  are  compared.  In  §3.1  the 
evolution  of  the  turbulence  kinetic  energy  across  the  shock  are  discussed.  The  inflnen/v^ 
of  the  shock  on  the  vorticity  fluctuations  is  considered  in  §3.2.  This  is  followed  in 
§3.3  by  a  phjrsical  argument  for  the  role  of  the  upstream  entropy  fluctuations.  In  §3.4 
a  scaling  for  the  amount  of  vorticity  produced  by  the  upstream  entropy  fluctuations 
is  derived.  Then  in  §3.5  the  evolution  of  the  thermodynamic  fluctuations  and  the 
validity  of  Morkovin’s  hypothesis  across  the  shock  are  discussed. 

3.1.  Turbulence  kinetic  energy 

The  presence  of  upstteam  entropy  fluctuations  has  a  noticeable  effect  on  the  evolution 
of  turbxflence  kinetic  energy  across  the  shock  wave.  Figures  5  and  6  show  the 
^eanmse  evolution  of  turbulence  kinetic  energy  in  cases  129A  and  129B.  Note  that 

p/2  ^  ^2  behind  the  shock.  The  intermittency  associated  with  shock ’bscillation  is 
seen  to  cause  high  fluctuation  levels  in  the  vicinity  of  the  shock  (Debieve  &  Lacharme 
1986;  Lee  et  al.  1992).  The  width  of  this  intermittent  region  (denoted  by  nearly 
equals  the  ‘m^  shock  thickness’.  Using  the  mean  velocity  profile  to  determine  its 
value,  kifbjfuer  is  approximately  0.3  and  0.4  in  cases  1.29 A  and  1.29B  respectively^ 
We  focw  our  attention  on  the  evolution  of  kinetic  energy  outside  this  intermittent 
region  in  the  following  paragraphs. 

As  shown  in  figure  6,  kinetic  energy  levels  behind  the  shock  wave  are  noticeably 
higher  in  case  1.29B.  The  streamwise  component  is  affectedjaore  than  the  transverse 
components.  Comparison  of  the  peak  in  the  profile  of  a'^  (at  Jlcox  =  14)  reveals 
20  ^  higher  levels  in  case  1.29B.  Comparison  of  v'^  at  the  same  location  shows 

t  pie  smallest  njcsh  spacing  is  in  the  vicinity  of  the  shock,  and  is  equal  to  0.02/ko  in  both 
smdations.  the  shock-capturing  scheme  is  applied  over  a  width  equal  to  1.8/*#  on  either  ride 
Of  the  mesui  shock. 
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Figure  5.  Turbulence  kinetic  energy  from  DNS  of  the  Mach  1.29  shock  wave.  M  the  curves  are 

normalized  by  their  value  immediately  upstream  of  the  shock  {kox  =  8.97): - ,  case  1.29A; 

- ,  v'^,  case  1.29A;  • ,  case  1.29B;  *  ,  case  1.29B. 


Figure  6.  Evolution  of  turbulence  kinetic  energy  behind  a  Mach  1.29  shodc  wave  as  predicted  by 
linear  analysis.  All  the  curves  are  normalized  by  thqr  value  immediately  upstream  of  the  shock: 

- ,  without  entropy  fluctuations; - ,  without  entropy  fluctuatitnis;  •.*>ii'3,-v»ith 

entropy  fluctuations;  >  ,  with  entropy  fluctuations. 

’N-.  .  ' 

the  level  of  enhancement  to  be  a  modest  7%.  This  .enhancement  in  the  presence 
of  entropy  fluctoations  is  predicted  by  linear  analysis.  Note  that  quantitative  com¬ 
parison  of  turbulence  kinetic  energy  between  analysis  and  DNS  is  made  difficult 
by  the  viscous  decay  in  the  simulations.  If  the  peak  levels  of  kinetic  energy  be¬ 
hind  the  shock  wave  are  compared,  the  linear  an^sis  predictions  are  seen  to  be 
higher.  As  seen  from  figure  6,  the  amplification  of  as  predicted  by  linear  analysis 
is  1.47  and  1.95  in  the  absence  and  presence  of  entropy  fluctuations  (that  satisfy 
Moikovin’s  hypothesis)  respectively.  Corresponding  values  for  are  predicted  to 
be  1.22  and  1.32  respectively.  Thus  linear  analysis  predicts  that  the  upstream  en¬ 
tropy  fluctuations  would  enhance  the  amplification  of  u'^  by  33%,  and  that  of  by 
about  8%. 
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Figure  7.  Turbulence  kinetic  energy  from  DNS  of  the  Mach  1.8  shock  wav^oth  the  curves  are 
nonnalized  by  their  value  immediately  upstream  of  the  shodc.  •  ,  1^;  * , 


Figure  8.  Evolution  of  turbulence  kinetic  energy  behind  a  Madi  1.8  shock  wave  as  predicted  by 
linear  anal^is.  All  the  curves  are  normalized  by  their  value  immediately  upstream  of  the  shock: 

>  without  entropy  fluctuations;  - - ,  without  entropy  fluctuations;  •  ,  with 

entropy  fluctuations;  * ,  j/*,  with  entropy  fluctuations. 

Figures  7  and  8  show  the  evolution  of  turbulence  kinetic  energy  across  the  Mach 
1.8  shock  wave.  The  amplification  levels  are  noticeably  higher  than  those  across  the 
Mach  U9  shock.  Using  values  me^ured  at  kox^  13.4  (the  location  behind  the  shock 
where  is  maximum)  ^nd  are  seen  to  amplify  by  2.25  and  1.61  respectively. 
These  levels  are  lower  than  predicted  by  linear  analysis,  where  Morkovin^s  hypothesis 
is  assumed  to  exactly  hold  upstream  of  the  shock.  The  predicted  (far-field)  values 
of  kinetic  energy  amplification  are  3.44  and  2.2  respectively.  Linear  analysis  predicts 
that  the  upstream  entropy  fluctuations  would  enhance  the  amplification  of  and 
by  110%  and  40%  respectively  across  a  Mach  1.8  shock. 

Note  that  ufT'  is  negative  if  Morkovin’s  hypothesis  is  valid.  Linear  analysis 
shows  that  the  sign  of  the  correlation  between  the  upstream  velocity  and  entropy 
fields  is  of  considerable  importance.  If  the  upstream  correlation  between  i/  and 
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Rguse  9.  (a)  and  (b)  in  the  far  field  of  the  shock  wave  as  predicted  by  linear  analysis.  All 

ttte  curves  are  normalized  by  their  upstream  values.  — —  ,  pure  vorticity; - -  ,  Morkovin’s 

hypothesis  satisfied  upstream; . ,  A,  =  0.54,^,  =  0*; - ,  Ar  *=  0.54,^,  =  180®. 


T  is  negative,  then  the  upstream  entropy  fluctuations  are  seen  to  enhance  kinetic 
energy  amplification.  This  is  the  trend  seen  in  the  DNS.  However,  if  vlV  is  positive 
upstream  of  the  shock,  the  amplification  of  kinetic  energy  is  strongly  suppressed. 
Practical  situations  where  this  might  occur  include  strongly  cooled  boundary  layers, 
and  mixing  layers  where  the  mean  velocity  and  density  gradients  are  of  opposite  sign. 
This  influence  of  entropy  fluctuations  increases  with  increasing  Mach  number.  Figures 
9{a)  and  9(h)  show  linear  analysis  predictions  of  kinetic  energy  amplification  in  the 
far  field  of  the  shock  wave.  If  Morkovin’s  hypothesis  is  assumed  to  hold  upstream, 
the  amplification  of  u'^  is  predicted  to  be  enhanced  by  more  than  100%  for  Mach 
numbers  exceeding  about  1.8.  The  primary  reason  for  this  increase  in  enhancement 
with  increasing  Mach  number  is  that  iPnta/p)/{qTm/U)  increases  with  Mach  number 
if  Morkovin’s  hypothesis  is  satisfied  upstream  of  the  shock  wave. 

3.2.  Vorticity  fluctuations 

All  components  of  vorticity  are  affected  by  the  upstream  entropy  fluctuations.  As 
shown  in  figure  10(a),  the  levels  of  vorticity  fluctuations  behind  the  shock  wave  are 
higher  in  case  1.29B.  The  amplification  of  <o^  is  seen  to  increase  by  8.7%  when 
entropy  fluctuations  are  present  upstream  of  the  shock.  This  increase  is  qualitatively 
predicted  by  linear  analysis.  The  increase  in  amplification  predicted  by  analysis  is 
much  higher  -  about  19.4%.  The  primary  reason  for  this  difference  betweeno^alysis... 
and  DNS  is  believed  to  be  the  strict  imposition  of  Morkovin’s  hypothesis  upstream 
of  the  shock  in  the  analysis.  Owing  to  the  absence  of  mean  temperature  gradients, 
the  upstream  fluctuations  in  the  simulation  only  approximately  satisfy  Morkovin’s 
hypothesis  (table  1).  Support  for  this  reasoning  is  provided  by  the  fact  that  the 
vorticity  amplification  in  case  l^A  is  within  6.3%  of  analysis  while  the  deviation  in 
case  1.29B  is  about  16.2%. 

Linear  analysis  predicts  the  transverse  vorticity  components  to  arnplify  across  the 
shock  wave,  and  remain  constant  downstream.  The  streamwise  vorticity  component 
is  predicted  to  remain  constant  across  and  downstream  of  the  shock.  DNS  shows  that 
the  streamwise  vorticity  component  does  indeed  remain  constant  across  the  shock. 
However,  it  increases  immediately  downstream  of  the  shock  wave.  In  fact,  the  peak 

level  of  behind  the  shock  wave  in  case  129B  is  about  7.1%  greater  than  its 
upstream  value.  Lee  et  al.  (1992)  noted  for  vortical  upstream  turbulence  that  this 
behaviour  is  nonlinear  in  nature;  it  is  caused  by  the  stretching  of  vorticity  fluctuations 
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Figure  10.  Streamwise  evolution  of  vorticity  fluctuations  from  DNS  of  (a)  the  Marh  129  shock 
wave  and  (h)  the  Mach  1.8  shock  wave.  All  the  curves  are  normaliad  by  their  value  immediately 

upstream  of  th^ock  (ktpc  =  8.97).  (^|- - case  1.29A; - ,  case  1.29A;  •  . 

case  U9B;  * ,  case  U9B.  (b)  • ,  d^^■,  k  , 


by  the  fluctuating  strain  rate.  In  the  presence  of  upstream  entropy  fluctuations,  it 
is  likely  that  nonlinear  baroclinic  effects  might  provide  an  additional  contribution. 
This  downstream  rise  in  is  noticeably  higher  behind  the  Mach  1.8  shock.  As 

shown  in  figure  10(6),  the  peak  level  of  behind  the  Mach  1.8  shock  is  about  twice 
its  upstream  value.  By  comparison,  Lee  et  al-'s  results  (in  the  absence  of  upstream 

entropy  fluctuations)  show  that  behind  a  Mach  2  shock  is  about  1.5  times  its 

upstream  value.  The  inability  to  predict  this  rise  in  (o'j^  behind  the  shock  wave  is 
probably  the  most  serious  limitation  of  the  linear  analysis.  Also,  it  is  interesting  to 
note  (figures  7  and  10(6))  that  the  vorticity  fluctuations  in  the  DNS  approach  isotropy 
more  rapidly  than  the  velocity  fluctuations.  This  is  consistent  with  the  notion  that 
the  small  scales  would  return  to  isotropy  more  rapidly  than  the  larger  scales.  Linear 
analysis  is  used  to  determine  the  influence  of  mean  Mach  number  on  the  amplification 
of  the  transverse  vorticity  fluctuations  in  figure  11.  As  observed  with  the  turbulence 
kinetic  energy,  negative  upstream  vlT  is  seen  to  enhance  the  amplification  of  vorticity 
fluctuations,  while  positive  correlation  is  seen  to  suppress  it 

3.3.  A  simple  explanation 

An  explanation  is  provided  here  for  the  influence  of  entropy  fluctuations  on  the 
evolution  of  a  turbulent  flow  across  a  shock  wave.  Note  that  HayesM1957)  results 
(equation  (30)  in  his  paper)  for  the  jump  in  vorticity  across  an  unsteady  shodc  wave 
require  the  unsteady  evolution  of  the  shock  front,  to  be  known  accurately  enough  that 
derivatives  normal  and  tangent  to  the  front  may  be  obtained.  This  is  quite  dfficult  to 
obtain  from  a  computation  where  the  shock  wave  is  ‘captured’  in  the  Navier-Stokes 
equations.  It  is  therefore  not  i^ssible  to  use  Hayes’  results  to  interpret  the  DNS. 
Consider  the  following  idealization  of  the  mean  field  associated  with  the  shock: 

£7  =  l/(x),  p  =  p(x),  p=p(x).  (3.1) 

Linear^g  the  Euler  equations  about  the  above  mean  flow  yields  the  following 
governing  equations  for  the  vorticity  fluctuations  (denoted  by  co'): 


J.  TTm'  — 


Py,  .  K  _ 


/ 


/ 
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Figure  1 1.  Amplification  of  across  the  shock  wave  as  predicted  by  linear  analysis. - , 

Pure  vorticity ; - ,  Morkovin’s  hypothesis  satisfied  upstream; . ,  A  =  0.54,  =  0“ ; - , 

A,  =0.54,  ^,  =  180". 

Note  that  although  irrotational  acoustic  fluctuations  are  generated  at  the  shock, 
their  contribution  to  the  downstream  kinetic  energy  is  not  significant;  linear  analysis 
(Mahesh  et  al.  1996)  shows  that  vortical  fluctuations  account  for  about  99%  of 
the  turbulence  kinetic  energy  behind  the  shock  wave.  Equation  (3.2)  represents  the 
effects  of  bulk  change  on  the  evolution  of  vorticity  fluctuations  across  the  shocL 
Linear  effects  due  to  shock  distortion  are  absent  Also,  the  effects  of  shock  curvature 
(Hayes  1957)  are  ignored.  As  a  result  of  these  assumptions,  the  relations  obtained  are 
expected  to  function  more  as  scaling  laws  than  predictive  formulae. 

If  the  incident  disturbance  comprises  vorticity  and  entropy  fluctuations,  then  p'  =  0 
in  the  incident  field  and  hence 

=  (33) 

The  term  —co'Ux  represents  the  effect  of  bulk  compression.  The  drop  in  mean  velocity 
across  the  shock  wave  indicates  that  this  term  would  enhance  the  incident  vorticity 
fluctuations.  The  incident  entropy  fluctuations  produce  vorticity  at  the  shock  wave 
through  die  baroclinic  term.  The  baroclinic  contribution  can  enhance  or  expose  the 
effect  of  bulk  compression.  The  phase  difference  between  the  upstream  vorticity  and 
entropy  waves  determines  whether  enhancement  or  opposition  is  observetL 

Consider  for  example  the  plane  vortidty-entropy'wave  represented  by  (A  Ic-e).  It 
is  easily  shown  that 

-co’Ux-^Px'-A,UUx-A,l^.  (3.4) 

r  p 

Since  t/,  is  negative  and  p^  is  positive  across  a  shock  wave,  the  two  sources  of 
vorticity  are  of  the  same  sign  if  At  and  A,  are  of  the  same  si^.  They  oppose  each 
other  if  A*  and  A,  are  of  opposite  sign.  Thus  if  vi  and  T'  are  negatively  correlated,  the 
entropy  field  enhances  the  amplification  of  fluctuating  vorticity.  On  the  other  hand,  a 
positive  correlation  between  vl  and  T  suppresses  the  amplification  of  vorticity  across 
the  shock  wave. 

Further  insight  is  gained  from  a  schematic  illustration  of  this  effect  Figure  12  shows 
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Bulk  cotni»ession 


Baroclinic 


Shock  wave 

Figure  11  A  spherical  clement  of  fluid  passing  through  a  shock  wave.  The  effects  of  bulk 
compression  and  baroclinic  vorticity  production  are  shown,  t/  and  T'  are  negatively  correlated 
upstream  of  the  shock  wave. 


a  fluid  element  of  circular  cross-section  passing  through  a  shock  wave.  The  geometric 
centre  of  this  element  is  denoted  by  Cf,  while  Cm  denotes  the  centre  of  mass.  The 
disturbance  field  associated  with  the  fluid  element  is  that  of  a  vorticity— entropy  wave. 
The  element  therefore  exhibits  solid-body  rotation  (with  associated  vorticity  o')  which 
is  assumed  positive  in  the  direction  shown.  Also,  the  density  gradient  associated  with 
the  entropy  wave  causes  the  centre  of  mass  to  differ  from  the  centre  of  force  (the 
geometric  centre).  Note  that  Cm  is  below  Cp  if  the  correlation  between  i/  and  T' 
is  negative.  Bulk  compression  compresses  the  element  in  the  streamwise  direction 
thereby  enhancing  the  rotation.  In  addition,  the  shock  wave  exerts  a  pressure  force 
(associated  with  the  adverse  pressure  gradient)  that  passes  through  Cp.  This  pressure 
force  would  exert  a  torque  about  the  centre  of  mass.  This  torque  manifests  itself  as 
the  baroclinic  source  of  vorticity.  Note  that  if  C*/  is  below  Cp  the  baroclinic  rotation 
is  in  the  s^e  direction  as  the  rotation  due  to  bulk  compression.  It  is  in  the  opposite 
direction  if  Cm  is  above  Cp  (positive  correlation  between  if  and  T').  The  upstream 
correlation  between  if  and  T'  thus  determines  the  location  of  Cm  with  respect  to 
Cp,  and  thereby  the  relative  sense  of  rotation  that  the  baroclinic  torque  produces. 


3.4.  Scaling  of  the  evolution  of  vorticity  across  a  shock  yvdve 

Equation  (3.3)  is  used  to  derive  approximate  expresaons  for  the  evolution  of  vorticity 
fluctuations  across  the  shock.  The  expressions  are  evaluated  by  comparing  to  the 
linear  analysis  predictions.  Equation  (3  J)  is  rewritten  as 

DtCt;©')  = -0  t/Px  (3.5) 

where  D,  denotes  the  material  derivative  5/flr  +  Vd/dx.  Using  the  relation  o-  = 
^pUUx  in  the  above  equation  yields 
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Figure  13.  (a)  Ifascaled  and  {Ja)  scaled  rjn.s.  vorticity  produced  in  the  interaction  of  an  entropy 
wave  with  a  shock  wave. - (Mi  =  1.25), - (Mi  =  1 J), . (Mi  =2), - (Mi  *  2.5). 


The  shock  wave  is  approximated  as  a  discontinuity.  {U\  is  expressed  as  A(U^)5(x), 
where  A(17^)  repre^ts  the  dilSerence  in  across  the  shock  wave,  and  5(x)  denotes 
the  Dirac  delta  function.  An  approximate  solution  to  the  above  equation  is  obtained 
by  setting  p'y  equal  to  its  upstream  value  (see  (A  Ic-e)).  Transforming  coordinates  to 
x!  —  x—Vt,x  ^  t  and  integrating  yields  the  following  expression  for  the  change  in 
vorticity  across  the  shock  wave: 

Uzta^  -  UiO)[  ~  (3.7) 

which  yields 

Ucl  1  — “ 

(0'2^r<o\  +  —  AeUi^—f-  (3.8) 

where  r  =  U1/U2  is  obtained  from  the  Rankine-Hugoniot  equations.  Equation  (3.8) 
suggests  that  incident  vorticity  fluctuations  amplify  by  an  amount  equal  to  the  mean 
density  ratio.  The  vorticity  produced  by  the  incident  entropy  fluctuations  is  predicted 
to  scale  as  fc^sinvi(l  — 

These  expressions  are  next  compared  to  linear  ansdysis.  Figures  13(a)  and  13(fe) 
show  the  levels  of  vorticity  produced  when  an  entropy  wave  interacts  with  a  shocl^ 
Mean  Mach  numbers  from  1.25  to  2.5  are  considered.  Both  unsealed  and  scaled 
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FicmE  14.  The  intensities  of  thermodynamic  fluctuations  (a)  from  DNS  and  (b)  as  predicted 
by  linear  analysis,  for  case  1.29B.  All  variables  are  non-  dimensionalized  by  the  value  of  Prm/p 

immediately  upstream  of  the  shock  wave  (kox  =  8.97). -  -  T  /T  . 

(7-l)A/2a™/ir; - ,p^/yp.  ’ 


values  of  rjn.s.  vorticity  behind  the  shock  wave  are  plotted  as  a  function  of  incidence 
angle.  The  proposed  scaling  seems  to  yield  reasonable  collapse  of  the  curves  in  the 
propagating  regime.  The  validity  of  the  scaling  for  a  turbulent  field  would  depend 
upon  the  fraction  of  incident  waves  in  the  propagating  regime. 

The  interaction  of  a  vorticity  wave  with  a  shock  has  been  considered,  and  the 
scahng  compared  to  linear  analysis  (Mahesh  et  al.  1996).  The  scaling  was  less 
satisfactory,  unsealed  values  ranging  from  1  to  6  were  collapsed  to  vary  from  0.5  to 
1.4.  The  ^plification  of  incident  vortidty  was  very  nearly  equal  to  the  tTipan  density 
ratio  for  iricidence  angles  near  zero;  however,  a  systematic  deviation  was  seen  with 
increa^g  incidence  an^e.  This  evolution  of  the  vorticity  waves  is  in  agreement  with 
Jacquin,  Cambon  &  Blin  (1993)  who  observe  large  discrepancies  at  the  higher  Mach 
numbers  when  the  effect  of  shock  distortion  is  not  represented  in  the  linear  limit 

3.5.  Thermodynamic  fluctuations  and  Morkomn‘s  hypothesis 
The  thermodynamic  fluctuations  behind  weak  shock  waves  were  noted  by  Lee,  Lele 
&  Moin  (1994)  to  be  nearly  isentropic.  The  thermo^amic  field  in  case  1.29A 
follows  this  trend,  Le.  Pms/yp,  Prmstp  and  Tr^l{;y  — 1)7  are  nearly  equal  over  the 
en^  domain;  However,  upstream  entropy  fluctuations  were  not  present  in  Lee  et 
aVs  computations.  As  might  be  expected  the  downstream  thermodyntfinic  field  is 
not  isentropic  when  upstream  entropy  fluctuations  are  present  Hgure  14(a)  shows 
the  streamwise  evolution  of  the  pressure,  density;  and  temperature  fluctuations  in 
^  129B.  The  quantity  (y  -  l)M^Ur,„/U  is  also  shown.  This  aUows  the  weak 
form  of  Morkovin’s  hypothesis  to  be  evaluated  across  the  shock.  The  corresponding 
predictions  made  by  linear  analysis  are  shown  in  figure  14(6). 

Go^  qualitative  agreement  is  observed  between  analysis  and  simulation.  The 
intensity  of  pressure  fluctuations  in  the  near  field  is  seen  to  be  comparable  to  that 
of  the  density  and  temperature.  However,  the  pressure  fluctuations  decay  behind  the 
shock  "wave,  causing  their  far-field  intensity  to  be  smaller.  Considerable  deviation 
from  Morkoi^’s  hypothesis  is  observed  in  the  near  field  behind  the  shock.  The  extent 
of  deviation  in  the  far  field  is  seen  to  be  smaller.  The  linear  analysis  predicts  a  larger 
deviation  from  the  hypothesis  in  the  far  field  than  the  DNS.  Owing  to  viscous  decay, 
the  notion  of  far  field  is  not  as  precise  in  the  DNS  as  it  is  in  the  analysis.  Linear 
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Figure  15.  Thermodynamic  fluctuations  (a)  from  DNS  and  (h)  as  predicted  by  linear  analysis,  for 
case  1.8.  All  variables  are  non-dimensionalized  by  the  value  of  pm/p  immediately  upstream  of  the 
shock  wave. - ,  prms/p; - ,  Tr^/T, .  ,  (y  -  l)Af^iw/C;; - ,  pnw/yp- 


analysis  suggests  (figure  14(fe))  that  the  far-field  values  are  attained  at  approximately 
fcox  =  17.  If  the  weak  form  of  Morkovin’s  hypothesis  is  evaluated  at  this  location  in 
the  computation,  it  shows  behaviour  comparable  to  analysis.  Interestingly  however, 
the  validity  of  the  hypothesis  in  the  DNS  is  seen  to  increase  with  distance  downstream 
of  this  location.  The  exact  cause  of  this  trend  is  not  known.  Also,  although  the  r.m.s 
values  of  the  terms  in  the  hypothesis  approach  each  other  behind  the  shock  in  the 
DNS,  t^T!/ihmsTrms  does  not  approach  —1.  It  drops  in  magnitude  across  the  shock, 
and  decreases  further  in  magnitude  downstream,  e.g.  the  correlation  coefficient  at 
fcox  =  20  is -0.54. 

A  similar  trend  is  observed  in  the  Mach  1.8  computation.  Morkovin’s  hypothesis 
is  evaluated  across  the  Mach  1.8  shock  wave  in  figures  15(a)  and  15(b)  respectively. 
The  far-field  values  in  the  analysis  are  attained  at  about  kox  —  15.  Both  analysis 
and  DNS  show  considerable  deviation  from  the  hypothesis  at  this  location.  However, 
the  intensities  of  the  density  and  temperature  fluctuations  decay  behind  the  shock 
wave  in  the  DNS,  resulting  in  increas^  validity  of  the  hypothesis  behind  the  shock. 
Interestingly,  although  priw/p,  T^/T  and  (y  —  l)M^Urm/U  approach  each  other 
behind  the  shock,  the  intensity  of  the  pressure  fluctuations  is  not  negligible.  As  seen 
in  figure  15(a),  prms/yp  is  nearly  of  the  same  magnitude  as  T„w/T  at  about  ^  40. 
Also,  ufT'/urasTrms  does  uot  approach  —1.  Its  value  at  /cqx  —  20  is  -^22, 

Linear  analysis  is  used  to  examine  the  influence  of  mean  Mach  number  on  the 
validity  of  Morkovin's  hypothesis  across  die  shock  ifrave.  The  incident  fluctuations 
are  constrained  to  satisfy  the  strict  form  of  the  hypothesis.  The  fluctuations  in  the  far 
field  are  then  exanoined  to  see  if  the  hypothesis  holds  in  the  rjn.s.  sense.  The  results 
(figure  16)  show  that  the  first  part  of  the  hypothesis,  Le.  prm/p  =  Tm^/T,  is  still  a 
good  approximation  behind  the  shock  wave  (especially  if  Afi  is  less  than  about  2). 
However,  the  part  of  the  hypothesis  that  relates  T'  to  t/  exhibits  large  deviation  with 
Mach  number.  This  behaviourjs  explained  below. 

The  equation  p'/P  =  — T'/T  is  obtained  by  setting  to  zero  in  the  linearized 
equation  of  state.  It  amounts  to  neglecting  the  acoustic  mode  in  comparison  to 
the  entropy  mode.  As  seen  from  {Aid),  it  holds  exactly  for  the  upstream  turbu¬ 
lence.  Upon  interaction  with  the  shock,  the  incident  fields  of  vorticity  and  entropy 
fluctuations  generate  acoustic  waves.  TTie  generation  of  acoustic  waves  is  however 
accompanied  by  amplification  of  the  incident  entropy  fluctuations.  Also,  a  fraction  of 
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Figure  16.  Evaluation  of  Morkt^’s  hypothesis  in  the  far  field  of  the  shock  wave  using  linear 

analysis.  AD  the  curves  are  non-dimensionalized  by  the  upstream  value  of  Ur—IJJ  _  o  tn- 

- .WT; . 


the  acoustic  waves  decays  behind  the  shock.  As  a  result,  the  acoustic  contribution  to 
the  thermod3maimc  fluctuations  in  the  far  field  is  significant  only  at  the  larger  Mach 
numbers.  The  firet  part  of  the  weak  form  of  Morkovin’s  hypothesis  is  therefore  a 
good  approximation  behind  shock  waves  of  moderate  strength. 

The  relation 

Y  =  -(y  -  ■  (3.9) 

is  obtained  by  assuming  that  stagnation  temperature  fluctuations  are  ■tman  in  the 
linear  limit.  Denote  the  stagnation  temperature  by  To, 


To  =  T  +  T'4 


2a 


(3.10) 


where  C,  denotes  the  specific  heat  at  constant  pressure.  Linearizing  the  above 
equation  yields  the  following  expression  for  fluctuations  in  stagnation  temperature: 
To  =  T  +  Ui/jCp.  Setting  Tq  to  zero  and  rearranging  yields  (3.9). 

In  the  linear  liimt,  fluctuations  in  stagnation  temperature  obey  the  relation, 
DTo/Dt  =  dp^/dt  ^hind  the  shock  wave  (D/Dt  denotes  the  material  derivative  based 
on  the  mean  velocity).  Decomposing  the  temperature  field  into  acoustic  aiid  entropic 
components,  and  the  velocity  into  acoustic  and  vortical  componerits  allovi^  Hecompo- 
ahon  of  T^into  T^  and  such  that  DT^^/Dr  =  0  and  DT^^/Dr  »  d^/dt. 
XMe^  the  Rankme-Hugomot  equations  may  be  used  to  show  that  (3.9)  cannot  be 
valid  behind  a  shock  wave  if  it  is  assumed  to  hold  upstream  of  the  shock.  The  energy 
equation  requires  the  stagnation  temperature  to  be  constant  across  the  sho<±  in  a 
frame  of  reference  that  moves  at  the  instantaneous  speed  of  the  shock  wave,  Le. 

+  ,3.„, 

p  2Cp 

Linearizing  the  above  equation  and  constraining  the  incident  fluctuations  to  satisfy 
(3.9)  yields 


Uz), 


(3.12) 


i.e.  the  fluctuations  of  stagnation  temperature  are  not  zero  immediately  behind  the 


Figure  17.  Decomposition  of  stagnation  temperature  fluctuations  using  linear  analysis.  All  the 

curves  are  non-dimensionalized  by  the  upstream  value  of  Ur^/U. - ,  Far-field  value; - , 

contribution  due  to  vorticity  and  entropy  fluctuations; . ,  near-field  value. 


shock  wave.  Dividing  through  by  Tz  and  rearranging,  we  get 

5+(r-l)M||  =  -(,-lW=(§-l)|.  (3.13) 

Thus,  applicability  of  Morkovin’s  hypothesis  immediately  behind  the  shock  wave 
requires  that 

(r-i)M,(^-i)i  ~  0.  (3.14) 

The  rms.  values  of  the  near-field  {kox  =  10)  stagnation  temperature  are  plotted 
in  figure  17  as  a  function  of  the  mean  Mach  number.  The  plotted  values  are 
seen  to  be  comparable  to  the  terms  in  Morkovin’s  hypothesis  (figure  16).  The 
hypothesis  is  therefore  invalid  immediately  behind  the  shock  wave.  The  above 
argument  is  nest  extended  to  show  why  the  hypothesis  does  not  hold  in  the  far  field. 
Decomposition  of  the  stagnation  temperature  fluctuations  into  vorticity-entropy  and 
acoustic  components  shows  that  both  near-  and  far-field  values  of  the  stagnation 
temperature  fluctuations  are  dominated  by  the  vorticity-entropy  component  As 
shown  in  figure  17,  the  vorticity-entropy  component  is  nearly  equ^  to  the  total  level 
in  the  far  field,  while  its  contribution  to  the  near-field  level  is  greater  than  80% 
over  the  range  of  Mach  numbers  shown.  Figure  17  and  (3.13)  therefore'^skow  that 
appreciable  levels  of  stagnation  temperature  fluctuations  are  generated  immediately 
behind  the  shock  wave  due  to  oscillation  of  the  shodC'front  Most  of  these  fluctuations 
arise  from  vorticity-entropy  fluctuations,  which  convect  downstream  to  generate  an 
appreciable  level  of  stagnation  temperature  fluctuations  in  the  far  field  of  the  shock 
wave.  This  leads  to  inapplicability  of  Morkoyin’s  hypothesis  in  the  far  field. 


4.  Summaiy. 

Direct  numerical  shnulation  and  invisdd  linear  analysis  were  used  to  study  the 
interaction  of  a  normal  diock  wave  with  an  isotropic  turbulent  field  of  vorticity 
and  entropy  fluctuations.  Shock  waves  of  strength  Mach  129  and  Mach  1.8  w^ 
computed  using  DNS,  while  the  linear  analysis  consideied  a  range  of  Mach  numbers 
from  1  to  3.  Our  objective  was  to  study  the  role  of  the  upstream  entropy  fluctuations 
in  the  evolution  of  the  turbulent  flow  across  the  shock. 
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The  upstreM  entropy  flucmations  were  shown  to  significantly  affect  shock/turbu¬ 
lence  interaction.  The  magnitude  of  the  entropy  fluctuations,  and  the  sign  of  the 
upstream  yelodty-temperature  correlation  were  both  seen  to  be  important.  Higher 
levels  of  kinetic  energy  and  vorticity  amplification  were  observed  across  the  shock 
when  tt'  and  T'  were  negatively  correlated  upstream  of  the  shock.  Positive  correlation 
had  the  opposite  effect.  An  explanation  was  provided  to  explain  these  trends.  The 
evolution  of  fluctuating  vorticity  across  the  shock  wave  was  noted  to  have  two  impor¬ 
tant  contributions:  bulk  compression  of  incident  vorticity  and  baroclinic 'production 
of  vortidty  ^ough  the  incident  entropy  fluctuations.  The  upstream  correlation  be¬ 
tween  vorticity  and  entropy  fluctuations  was  shown  to  determine  whether  these  two 
sources  of  vorticity  enhance  or  oppose  each  other,  thereby  determining  kinetic  energy 
levels  behind  the  shock  wave.  A  scaling  was  then  proposed  for  the  evolution  of 
vorticity  across  the  shock  wave.  Since  Morkovin’s  hypothesis  is  known  to  apply  in 
adiabatic  turbulent  boundary  layers,  the  results  suggest  that  the  entropy  fluctuations 
in  the  boundary  layer  would  play  a  very  significant  role  in  the  interaction  between 
the  boundary  layer  and  a  shock  wave. 

T^e  v^dity  of  Morkovin’s  hypothesis  behind  a  shock  was  examiiiftfl  Linear 
analysis  indicates  that  neglecting  the  acoustic  mode  is  a  good  approximation  in 
the  far  field  of  shock  waves  of  moderate  strength  (Mi  <  2).  The  part  of  the 
hpjthesis  relating  t/  and  T'  was  predicted  to  be  invalid.  Non-negligible  oscillation 
of  the  shock  front  was  shown  to  be  responsible.  The  thermodynamic  fluctuations  in 
the  simulations  followed  linear  analysis  predictions  immediately  behind  the  shock. 
Interestingly  however,  the  terms  in  the  weak  form  of  the  hypothesis  approached  each 
other  as  the  solution  decayed  behind  the  shock.  Despite  this  trend,  the  pressure 
fluctuations  in  the  DNS  were  not  negligible  in  the  far  field.  Also,  IFF  did  not 
approach  -1.  Its  value  in  the  far  field  was  about  -0.54  and  -022  for  the  Mach  1.29 
and  Mach  1.8  shock  waves. 

This  study  was  supported  by  the  Air  Force  Office  of  Scientific  Research  under 
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twhnical  monitor.  The  authors  would  also  like  to  express  their  gratitude  to  NAS  and 
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Appendix  - 

The  linear  analysis  is  summarized  below.  The  interaction  of  a  sin^  vorticity- 
enteopy  wave  with  a  shock  is  first  considered  in  §A>1.  The  analysis  is  then  extended 
m  §A2  to  describe  the  evolution  of  the  turbulent  field  across  the  shock. 

Al.  Interaction  of  a  shock  with  a  plane  vorticity— entropy  wave 

The  two-dimensional  interaction  of  a  shock  wave  with  a  plane  vorticity-entropy  wave 
IS  schematicaUy  illustrated  in  figure  18.  Note  that  the  shock  wave  is  stationary  in 
the  mean.  The  variables  V,p,p,  T  and  M  denote  the  mean  velocity,  pressure,  density, 
temperature  and  Mach  number  respectively  and  subscripts  1  and  2  denote  the 
ups^am  and  downstream  states.  The  flow  upstream  of  the  shock  wave  is  perturbed 
by  me  weak  disturbance  field  of  the  incident  vortidty-entropy  wave  which  is  assumtHi 
tojje  a  plane  wave  that  makes  angle  ipi  with  the  x-axis.  The  variables 
and  T  represent  the  fluctuating  velocities,  pressure,  density  and  temperature.  ITie 
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Figure  18.  Schematic  of  the  interaction  of  a  vorticity-entropy  wave  with  a  shock  wave. 


incident  jSeld  has  the  following  form: 

^  =  lA,  ^  =  -mA,  (A  laj}) 

Ux 

p\=0,  (Alc-e) 

Pi  Tt  Pi 

where  m  =  cosv’i  and  I  =  sin  The  shock  wave  deforms  in  response  to  the  incident 
disturbance;  the  displacement  is  denoted  by  x  = 

Details  of  the  analysis  may  be  found  in  Mahesh  et  oL  (1996).  .  The  solution 
obtained  is  summarized  below.  The  solution  behind  the  shock  wave  is  a  superposition 
of  vortical,  acoustic  and  ^tropic  components.  It  has  the  following  dimensionless 
form: 


_  P  gpcXgMly-mVit)  1  Mnrx+ly-mUit) 

A,Ui  ^ 

—jj  ^mrx+fy-mU,t)^ 

Ajf  p2 

■  .  —8l  —  ^  j.  Q  ^iirx+ly-mUit) 

A„P2  7 

^  _  y  ^tkx^ly-mVit}  _Q  ^mx^y—mUi^ 

AfTz  7 


(A2fl) 

(A2h) 

(A2c) 


(A2d) 


(A2c) 


The  boundary  conditions  across  the  shock  wave  yield  the  following  expressions  for 
the  velocity  and  slope  of  the  shock  front: 

^  1  Ap  ftl  , 
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Figure  19.  Coordinate  system  used  in  the  interaction  of  isotropic  turbulence  with  a  shock  wave. 


The  variable  r  denotes  Ae  mean  density  ratio  pj/p,  across  the  shock  wave;  k(MutPt) 
represents  the  streamwise  wavenumber  of  the  acoustic  component  behind  the  shock 
wave.  It  IS  real_^m  the  propagating  regime  and  complex  in  the  decaying  regime.  The 
coefficients  F,H  and  K  are  associated  with  the  acoustic  component  The  vortical 
component  b  reprinted  by  G  and  7  while  Q  represents  the  entropy  component  The 
coeffidents  K,G,  I  and  Q  are  functions  of  Mi,  ipi  and  the  amplitude  ratio  (Ar) 
Md  ph^e  difference  {<f>r)  between  the  vorticity  and  entropy  waves,  Le.  y4,e^  =  A,/ A 
Expressions  for  the  coefficients  are  given  in  Mahesh  et  al.  (1996).  ' 


A,2.  Interaction  of  a  shodc  wave  with  an  isotropic  turbulent  field 

As  suggest^  by  Morkovin’s  hypothesis,  the  upstream  turbulence  is  represented  as 
a  random  three^ensional  field  of  vorticity-entropy  waves.  The  turbulent  velocity 

eld  IS  ^sumed  to  be  isotropic.  The  incident  turbulent  field  is  represented  as  a 
mperposition  of  plane  vorticity-entropy  waves  (Fourier  modes)  in  three  dimensions 
Each  of  these  wavw  would  mteract  independentiy  with  the  shock  wave  under  linear 
analysis.  For  a  given  upstream  spectrum,  the  interaction  of  each  wave  with  the 
shock  wave  is  predicted.  The  solution  is  integrated  over  aU  incident  wavM  to  obtain 
turbulence  statistics  behind  the  shock. 

The  three-dimensional  problem  k  related  to  the  two-dimensional  analysis  of  the 
precalmg  se(^on  as  follows.  Consider  an  incident  plane  wave  in  three  dimensions. 
As  shown  m  figure  19,  the  wavenumber  vector  of  the  wave  lies  in  a  plane  that  makes 
an^e  with  the  y-axis.  In  this  plane,  which  we  caff  the  (x,x,)-plane.  the  wave 
^es  angle  vi  with  the  x-axis.  It  is  readily  seen  that  the  (x,x,)-plane  is  identical 
to  ffie  plMe  of  mterachon  m  the  two-dimensional  problem.  The  solenoidal  nature 
of  ffie  madent  veloaty  field  requires  the  velocity  vector  of  file  wave  to  be  normal 
o  me  wavenumber  vector.  The  velocity  field  may  therefore  be  expressed  as  a  sum 
ot  two  components:  one  normal  to  the  wavenumber  vector  in  the  (x,Xr>pIane  and 
me  other  normd  to  the  (x,Xr)-plane  (the  (direction).  It  is  intuitively  clear  that  the 
^-component  of  velocity  would  pass  unchanged  through  the  shock  wave.  As  a  result, 
the  three-dimension^  problem  may  be  solved  using  results  of  the  two-dimensional 
analysis  m  the  (x,Xr)-plane. 
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Results  from  §A.l  are  used  to  obtain  expressions  for  the  energy  spectra  behind 
the  shock  wave.  The  spectra  depend  upon  the  upstream  three-dimensional  energy 
spectrum,  E{k).  This  paper  assumes  the  following  form  for  E{k): 


E{k) 


~(0 


y--2{klhf 


(A4) 


The  results  also  depend  upon  the  upstream  density  spectrum.  The  quantity  At/A„ 
may  be  represented  as  A,^',  where  At  and  <j>r  are  both  functions  of  the  wavenum¬ 
ber  vector.  Appropriate  functional  dependencies  may  be  assumed  depending  upon 
the  flow  being  considered.  This  paper  presents  resiflts  for  <l>r  —  0  and  it,  Le.  the 
density  field  is  either  perfectly  correlated  or  perfectly  anti-coirelated  with  the  ve¬ 
locity  field.  Also,  two  forms  of  the  upstream  density  spectrum  are  considered.  One 
case  assiunes  the  density  field  to  be  isotropic  with  the  same  three-dimensional  spec¬ 
trum  as  the  velocity  field,  ie.  Ar  is  assumed  constant  For  this  case  it  is  easily 
shown  that 


A 


(A5) 


The  second  case  assumes  that  the  density  field  satisfies  Morkovin’s  hypothesis  at 
every  wavenumber  If  the  velocity  fluctuations  are  isotropic,  it  is  easily  seen  that  the 
resulting  density  field  is  axisymmetiic,  Le. 


Ar  =  (7  —  l)M^sinv?i. 


(A6) 


The  spectra  are  then  numerically  integrated  to  obtain  the  turbulence  statistics  behind 
the  shock.  Further  details  are  provided  by  Mahesh  et  aL  (1996). 
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